问题：

C语言中三角函数的单精度变元约简

段干宜

2023-03-14

我已经实现了三角函数(sin，cos，arctan)的一些近似值，用单精度（32位浮点）计算，它们的精度约为+/-2ULP。

我的目标设备不支持任何或方法。它不提供FMA，而是提供MAC ALU。ALU和LU以32位格式计算。

我现在的目标是为正弦和余弦提供更大的输入范围（尽可能大，理想情况下是[FLT_MIN，flt_max])，这导致我减少参数。

我目前正在阅读不同的论文，比如K.C.Ng写的关于大参数的参数缩减：Good to the Last Bit，或者是关于这个新的参数缩减算法的论文，但我无法从中导出一个实现。

另外，我想提到两个stackoverflow问题，它们涉及到相关问题：有一个使用matlab和C++的方法，它基于我链接的第一篇论文。它实际上使用了matlab，cmath方法，并将输入限制在[0,20.000]。另一个已经在评论中提到了。它是一种在C中实现sin和cos的方法，使用了我无法使用的各种C库。由于这两个帖子都已经有几年了，可能会有一些新的发现。

所以我的问题是：有没有另一种有效的方法来减少正弦和余弦的争论，获得单一精度，避免大的LUT？上面提到的论文实际上侧重于双精度，使用高达1000位数，这不适合我的用法。

实际上，我还没有在C中找到任何实现，也没有针对单一精度计算的实现，我将非常感谢任何类型的提示/链接/示例...

易弘亮

2023-03-14

下面的代码是基于前面的一个答案，在这个答案中，我演示了如何通过对大小较小的参数使用拆分常数的Cody-Waite方法，以及对大小较大的参数使用Payne-Hanek方法，对三角函数执行相当精确的参数缩减。关于Payne-Hanek算法的详细信息，请参见这里，关于Cody-Waite算法的详细信息，请参见我之前的答案。

这里我做了必要的调整，以适应Asker平台的限制，因为不支持64位类型，不支持融合的乘法-加法，并且Math.h中的helper函数不可用。我假设float映射到IEEE-754binary32格式，并且有一种方法可以将这样的32位float重新解释为32位无符号整数，反之亦然。我已经通过标准可移植习惯用法实现了这种重新解释，即通过使用memcpy()，但是对于未指定的目标平台，也可以选择其他合适的方法，例如内联汇编、特定于机器的内部特性或volatile联合。

由于这段代码基本上是我以前的代码到一个更严格的环境的一个端口，它可能缺乏专门针对该环境的从头设计的优雅。基本上，我已经用一些位旋转替换了Math.h中的frexp()helper函数，用32位整数对模拟了64位整数计算，用32位定点计算替换了双精度计算（比我预期的要好得多），并用未使用的等效函数替换了所有FMAs。

重新工作的Cody-Waite部分的论点减少需要相当多的工作。显然，在没有FMA可用的情况下，我们需要确保常数π/2的组成部分中有足够数量的尾随零位（最低有效位除外），以确保乘积是准确的。我花了几个小时来尝试找出一个特定的分裂，它可以提供准确的结果，但也可以将转换点推到佩恩-哈内克方法的尽可能高的位置。

当指定use_fma=1时，用高质量的数学库编译测试应用程序的输出应该类似于以下内容：

Testing sinf ...  PASSED. max ulp err = 1.493253  diffsum = 337633490
Testing cosf ...  PASSED. max ulp err = 1.495098  diffsum = 342020968

当use_fma=0时，精确度略有下降：

Testing sinf ...  PASSED. max ulp err = 1.498012  diffsum = 359702532
Testing cosf ...  PASSED. max ulp err = 1.504061  diffsum = 364682650

请注意，使用最严格的浮点设置和编译器提供的最高程度的IEEE-754来编译代码是很重要的。对于我用来开发和测试这段代码的Intel编译器，这可以通过使用/fp:strict编译来实现。此外，参考数学库的质量对于准确评估这种单精度代码的ulp误差至关重要。Intel编译器附带了一个数学库，它提供了双精度的基本数学函数，在HA（高精度）变体中，ulp误差略高于0.5ulp。使用多精度参考库可能更好，但会使我在这里慢得太多。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>   // for memcpy()
#include <math.h>     // for test purposes, and when PORTABLE=1 or USE_FMA=1

#define USE_FMA   (0) // use fmaf() calls for arithmetic
#define PORTABLE  (0) // allow helper functions from math.h
#define HAVE_U64  (0) // 64-bit integer type available
#define CW_STAGES (3) // number of stages in Cody-Waite reduction when USE_FMA=0

#if USE_FMA
#define SIN_RED_SWITCHOVER  (117435.992f)
#define COS_RED_SWITCHOVER  (71476.0625f)
#define MAX_DIFF            (1)
#else // USE_FMA
#if CW_STAGES == 2
#define SIN_RED_SWITCHOVER  (3.921875f)
#define COS_RED_SWITCHOVER  (3.921875f)
#elif CW_STAGES == 3
#define SIN_RED_SWITCHOVER  (201.15625f)
#define COS_RED_SWITCHOVER  (142.90625f)
#endif // CW_STAGES
#define MAX_DIFF            (2)
#endif // USE_FMA

/* re-interpret the bit pattern of an IEEE-754 float as a uint32 */
uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/* re-interpret the bit pattern of a uint32 as an IEEE-754 float */
float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/* Compute the upper 32 bits of the product of two unsigned 32-bit integers */
#if HAVE_U64
uint32_t umul32_hi (uint32_t a, uint32_t b)
{
    return (uint32_t)(((uint64_t)a * b) >> 32);
}
#else // HAVE_U64
/* Henry S. Warren, "Hacker's Delight, 2nd ed.", Addison-Wesley 2012. Fig. 8-2 */
uint32_t umul32_hi (uint32_t a, uint32_t b)
{
    uint16_t a_lo = (uint16_t)a;
    uint16_t a_hi = a >> 16;
    uint16_t b_lo = (uint16_t)b;
    uint16_t b_hi = b >> 16;
    uint32_t p0 = (uint32_t)a_lo * b_lo;
    uint32_t p1 = (uint32_t)a_lo * b_hi;
    uint32_t p2 = (uint32_t)a_hi * b_lo;
    uint32_t p3 = (uint32_t)a_hi * b_hi;
    uint32_t t = (p0 >> 16) + p1;
    return (t >> 16) + (((uint32_t)(uint16_t)t + p2) >> 16) + p3;
}
#endif // HAVE_U64

/* 190 bits of 2/PI for Payne-Hanek style argument reduction. */
const uint32_t two_over_pi_f [] = 
{
    0x28be60db,
    0x9391054a,
    0x7f09d5f4,
    0x7d4d3770,
    0x36d8a566,
    0x4f10e410
};

/* Reduce a trig function argument using the slow Payne-Hanek method */
float trig_red_slowpath_f (float a, int *quadrant)
{
    uint32_t ia, hi, mid, lo, tmp, i, l, h, plo, phi;
    int32_t e, q;
    float r;

#if PORTABLE
    ia = (uint32_t)(fabsf (frexpf (a, &e)) * 0x1.0p32f); // 4.29496730e+9
#else // PORTABLE
    ia = ((float_as_uint32 (a) & 0x007fffff) << 8) | 0x80000000;
    e = ((float_as_uint32 (a) >> 23) & 0xff) - 126;
#endif // PORTABLE
    
    /* compute product x * 2/pi in 2.62 fixed-point format */
    i = (uint32_t)e >> 5;
    e = (uint32_t)e & 31;

    hi  = i ? two_over_pi_f [i-1] : 0;
    mid = two_over_pi_f [i+0];
    lo  = two_over_pi_f [i+1];
    tmp = two_over_pi_f [i+2];
 
    if (e) {
        hi  = (hi  << e) | (mid >> (32 - e));
        mid = (mid << e) | (lo  >> (32 - e));
        lo  = (lo  << e) | (tmp >> (32 - e));
    }

    /* compute 64-bit product phi:plo */
    phi = 0;
    l = ia * lo;
    h = umul32_hi (ia, lo);
    plo = phi + l;
    phi = h + (plo < l);
    l = ia * mid;
    h = umul32_hi (ia, mid);
    plo = phi + l;
    phi = h + (plo < l);
    l = ia * hi;
    phi = phi + l;

    /* split fixed-point result into integer and fraction portions */
    q = phi >> 30;               // integral portion = quadrant<1:0>
    phi = phi & 0x3fffffff;      // fraction
    if (phi & 0x20000000) {      // fraction >= 0.5
        phi = phi - 0x40000000;  // fraction - 1.0
        q = q + 1;
    }

    /* compute remainder of x / (pi/2) */
#if USE_FMA
    float phif, plof, chif, clof, thif, tlof;
    phif = 0x1.0p27f * (float)(int32_t)(phi & 0xffffffe0);
    plof = (float)((plo >> 5) | (phi << (32-5)));
    thif = phif + plof;
    plof = (phif - thif) + plof;
    phif = thif;
    chif =  0x1.921fb6p-57f; // (1.5707963267948966 * 0x1.0p-57)_hi
    clof = -0x1.777a5cp-82f; // (1.5707963267948966 * 0x1.0p-57)_lo
    thif = phif * chif;
    tlof = fmaf (phif, chif, -thif);
    tlof = fmaf (phif, clof, tlof);
    tlof = fmaf (plof, chif, tlof);
    r = thif + tlof;
#else // USE_FMA
    /* record sign of fraction */
    uint32_t s = phi & 0x80000000;
    
    /* take absolute value of fraction */
    if ((int32_t)phi < 0) {
        phi = ~phi;
        plo = 0 - plo;
        phi += (plo == 0);
    }
    
    /* normalize fraction */
    e = 0;
    while ((int32_t)phi > 0) {
        phi = (phi << 1) | (plo >> 31);
        plo = plo << 1;
        e--;
    }
    
    /* multiply 32 high-order bits of fraction with pi/2 */
    phi = umul32_hi (phi, 0xc90fdaa2); // (uint32_t)rint(PI/2 * 2**31)
    
    /* normalize product */
    if ((int32_t)phi > 0) {
        phi = phi << 1;
        e--;
    }

    /* round and convert to floating point */
    uint32_t ri = s + ((e + 128) << 23) + (phi >> 8) + ((phi & 0xff) > 0x7e);
    r = uint32_as_float (ri);
#endif // USE_FMA
    if (a < 0.0f) {
        r = -r;
        q = -q;
    }

    *quadrant = q;
    return r;
}

/* 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 {
        /* Cody-Waite style reduction. W. J. Cody and W. Waite, "Software Manual
           for the Elementary Functions", Prentice-Hall 1980
        */
#if USE_FMA
        j = fmaf (a, 0x1.45f306p-1f, 0x1.8p+23f) - 0x1.8p+23f; // 6.36619747e-1, 1.25829120e+7
        r = fmaf (j, -0x1.921fb0p+00f, a); // -1.57079601e+00 // pio2_high
        r = fmaf (j, -0x1.5110b4p-22f, r); // -3.13916473e-07 // pio2_mid
        r = fmaf (j, -0x1.846988p-48f, r); // -5.39030253e-15 // pio2_low
#else // USE_FMA
        j = (a * 0x1.45f306p-1f + 0x1.8p+23f) - 0x1.8p+23f; // 6.36619747e-1, 1.25829120e+7
#if CW_STAGES == 2
        r = a - j * 0x1.921fb4p+0f;  // pio2_high
        r = r - j * 0x1.4442d2p-24f; // pio2_low
#elif CW_STAGES == 3
        r = a - j * 0x1.921f00p+00f; // 1.57078552e+00 // pio2_high
        r = r - j * 0x1.6a8880p-17f; // 1.08043314e-05 // pio2_mid
        r = r - j * 0x1.68c234p-39f; // 2.56334407e-12 // pio2_low
#endif // CW_STAGES
#endif // USE_FMA
        *q = (int)j;
    }
    return r;
}

/* Approximate sine on [-PI/4,+PI/4]. Maximum ulp error with USE_FMA = 0.64196
   Returns -0.0f for an argument of -0.0f
   Polynomial approximation based on T. Myklebust, "Computing accurate 
   Horner form approximations to special functions in finite precision
   arithmetic", http://arxiv.org/abs/1508.03211, retrieved on 8/29/2016
*/
float sinf_poly (float a, float s)
{
    float r, t;
#if USE_FMA
    r =              0x1.80a000p-19f;  //  2.86567956e-6
    r = fmaf (r, s, -0x1.a0690cp-13f); // -1.98559923e-4
    r = fmaf (r, s,  0x1.111182p-07f); //  8.33338592e-3
    r = fmaf (r, s, -0x1.555556p-03f); // -1.66666672e-1
    t = fmaf (a, s, 0.0f); // ensure -0 is passed through
    r = fmaf (r, t, a);
#else // USE_FMA
    r =         0x1.80a000p-19f; //  2.86567956e-6
    r = r * s - 0x1.a0690cp-13f; // -1.98559923e-4
    r = r * s + 0x1.111182p-07f; //  8.33338592e-3
    r = r * s - 0x1.555556p-03f; // -1.66666672e-1
    t = a * s + 0.0f; // ensure -0 is passed through
    r = r * t + a;
#endif // USE_FMA
    return r;
}

/* Approximate cosine on [-PI/4,+PI/4]. Maximum ulp error with USE_FMA = 0.87444 */
float cosf_poly (float s)
{
    float r;
#if USE_FMA
    r =              0x1.9a8000p-16f;  //  2.44677067e-5
    r = fmaf (r, s, -0x1.6c0efap-10f); // -1.38877297e-3
    r = fmaf (r, s,  0x1.555550p-05f); //  4.16666567e-2
    r = fmaf (r, s, -0x1.000000p-01f); // -5.00000000e-1
    r = fmaf (r, s,  0x1.000000p+00f); //  1.00000000e+0
#else // USE_FMA
    r =         0x1.9a8000p-16f; //  2.44677067e-5
    r = r * s - 0x1.6c0efap-10f; // -1.38877297e-3
    r = r * s + 0x1.555550p-05f; //  4.16666567e-2
    r = r * s - 0x1.000000p-01f; // -5.00000000e-1
    r = r * s + 0x1.000000p+00f; //  1.00000000e+0
#endif // USE_FMA
    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 with USE_FMA = 1: 1.495098  */
float my_sinf (float a)
{
    float r;
    int i;

    a = a * 0.0f + a; // inf -> NaN
    r = trig_red_f (a, SIN_RED_SWITCHOVER, &i);
    r = sinf_cosf_core (r, i);
    return r;
}

/* maximum ulp error with USE_FMA = 1: 1.493253 */
float my_cosf (float a)
{
    float r;
    int i;

    a = a * 0.0f + a; // inf -> NaN
    r = trig_red_f (a, COS_RED_SWITCHOVER, &i);
    r = sinf_cosf_core (r, i + 1);
    return r;
}

/* re-interpret bit pattern of an IEEE-754 double as a uint64 */
uint64_t double_as_uint64 (double a)
{
    uint64_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

double floatUlpErr (float res, double ref)
{
    uint64_t i, j, err, refi;
    int expoRef;
    
    /* ulp error cannot be computed if either operand is NaN, infinity, zero */
    if (isnan (res) || isnan (ref) || isinf (res) || isinf (ref) ||
        (res == 0.0f) || (ref == 0.0f)) {
        return 0.0;
    }
    /* Convert the float result to an "extended float". This is like a float
       with 56 instead of 24 effective mantissa bits.
    */
    i = ((uint64_t)float_as_uint32(res)) << 32;
    /* Convert the double reference to an "extended float". If the reference is
       >= 2^129, we need to clamp to the maximum "extended float". If reference
       is < 2^-126, we need to denormalize because of the float types's limited
       exponent range.
    */
    refi = double_as_uint64(ref);
    expoRef = (int)(((refi >> 52) & 0x7ff) - 1023);
    if (expoRef >= 129) {
        j = 0x7fffffffffffffffULL;
    } else if (expoRef < -126) {
        j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
        j = j >> (-(expoRef + 126));
    } else {
        j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
        j = j | ((uint64_t)(expoRef + 127) << 55);
    }
    j = j | (refi & 0x8000000000000000ULL);
    err = (i < j) ? (j - i) : (i - j);
    return err / 4294967296.0;
}

int main (void) 
{
    float arg, res, reff;
    uint32_t argi, resi, refi;
    int64_t diff, diffsum;
    double ref, ulp, maxulp;

    printf ("Testing sinf ...  ");
    diffsum = 0;
    maxulp = 0;
    argi = 0;
    do {
        arg = uint32_as_float (argi);
        res = my_sinf (arg);
        ref = sin ((double)arg);
        reff = (float)ref;
        resi = float_as_uint32 (res);
        refi = float_as_uint32 (reff);
        ulp = floatUlpErr (res, ref);
        if (ulp > maxulp) {
            maxulp = ulp;
        }
        diff = (resi > refi) ? (resi - refi) : (refi - resi);
        if (diff > MAX_DIFF) {
            printf ("\nerror @ %08x (% 15.8e): res=%08x (% 15.8e)  ref=%08x (%15.8e)\n", argi, arg, resi, res, refi, reff);
            return EXIT_FAILURE;
        }
        diffsum = diffsum + diff;
        argi++;
    } while (argi);
    printf ("PASSED. max ulp err = %.6f  diffsum = %lld\n", maxulp, diffsum);

    printf ("Testing cosf ...  ");
    diffsum = 0;
    maxulp = 0;
    argi = 0;
    do {
        arg = uint32_as_float (argi);
        res = my_cosf (arg);
        ref = cos ((double)arg);
        reff = (float)ref;
        resi = float_as_uint32 (res);
        refi = float_as_uint32 (reff);
        ulp = floatUlpErr (res, ref);
        if (ulp > maxulp) {
            maxulp = ulp;
        }
        diff = (resi > refi) ? (resi - refi) : (refi - resi);
        if (diff > MAX_DIFF) {
            printf ("\nerror @ %08x (% 15.8e): res=%08x (% 15.8e)  ref=%08x (%15.8e)\n", argi, arg, resi, res, refi, reff);
            return EXIT_FAILURE;
        }
        diffsum = diffsum + diff;
        argi++;
    } while (argi);
    diffsum = diffsum + diff;
    printf ("PASSED. max ulp err = %.6f  diffsum = %lld\n", maxulp, diffsum);
    return EXIT_SUCCESS;
}

C语言中三角函数的单精度变元约简

共有1个答案

相关问答

相关文章

相关阅读

相关工具

相关文档