c - 使用标准 C 数学库实现 sinpi() 和 cospi()

Question

函数 sinpi(x) 计算 sin(πx)，函数 cospi(x) 计算 cos(πx)，其中与 π 的乘法在内部是隐式的功能。这些函数最初由 Sun Microsystems 在 late 1980s 中作为扩展引入 C 标准数学库。 . IEEE Std 754™-2008 在第 9 节中指定了等效函数 sinPi 和 cosPi。

有许多自然发生 sin(πx) 和 cos(πx) 的计算。一个非常简单的例子是 Box-Muller 变换(G. E. P. Box 和 Mervin E. Muller，“关于随机正态偏差生成的注释”。数理统计年鉴，第 29 卷，第2, pp. 610 - 611)，给定两个均匀分布的独立随机变量 U₁ 和 U₂，生成标准正态分布的独立随机变量 Z₁ 和 Z₂:

Z₁ = √(-2 ln U₁) cos (2 π U₂)
Z₂ = √(-2 ln U₁) sin (2 π U₂)

另一个例子是度参数的正弦和余弦计算，如使用 Haversine 公式计算大圆距离:

/* This function computes the great-circle distance of two points on earth 
   using the Haversine formula, assuming spherical shape of the planet. A 
   well-known numerical issue with the formula is reduced accuracy in the 
   case of near antipodal points.

   lat1, lon1  latitude and longitude of first point, in degrees [-90,+90]
   lat2, lon2  latitude and longitude of second point, in degrees [-180,+180]
   radius      radius of the earth in user-defined units, e.g. 6378.2 km or 
               3963.2 miles

   returns:    distance of the two points, in the same units as radius

   Reference: http://en.wikipedia.org/wiki/Great-circle_distance
*/
double haversine (double lat1, double lon1, double lat2, double lon2, double radius)
{
    double dlat, dlon, c1, c2, d1, d2, a, c, t;

    c1 = cospi (lat1 / 180.0);
    c2 = cospi (lat2 / 180.0);
    dlat = lat2 - lat1;
    dlon = lon2 - lon1;
    d1 = sinpi (dlat / 360.0);
    d2 = sinpi (dlon / 360.0);
    t = d2 * d2 * c1 * c2;
    a = d1 * d1 + t;
    c = 2.0 * asin (fmin (1.0, sqrt (a)));
    return radius * c;
}

对于 C++，Boost 库提供了 sin_pi和 cos_pi , 并且一些供应商提供 sinpi 和 cospi 功能作为系统库的扩展。例如Apple在iOS 7中加入了__sinpi、__cospi以及对应的单精度版本__sinpif、__cospif， OS X 10.9(presentation，幻灯片 101)。但是对于许多其他平台，没有可供 C 程序轻松访问的实现。

与使用例如sin(M_PI * x)和cos(M_PI * x)，sinpi和cospi的使用提高了精度通过与 π 的内部 乘法减少舍入误差，并且由于更简单的参数减少还提供了性能优势。

如何使用标准 C 数学库以合理高效且符合标准的方式实现 sinpi() 和 cospi() 功能？

Answer 1

为简单起见，我将重点介绍 sincospi()，它同时提供正弦和余弦结果。然后可以将 sinpi 和 cospi 构造为丢弃不需要数据的包装函数。在许多应用程序中，浮点标志的处理(参见fenv.h)是不需要的，大多数时候我们也不需要errno错误报告，所以我将省略这些。

基本的算法结构很简单。由于非常大的参数总是偶数整数，因此是 2π 的倍数，因此它们的正弦和余弦值是众所周知的。在记录象限信息时，其他参数被折叠到 [-¼,+¼] 范围内。多项式 minimax approximations用于计算初级近似区间上的正弦和余弦。最后，象限数据通过结果的循环交换和符号变化将初步结果映射到最终结果。

正确处理特殊操作数(特别是 -0、无穷大和 NaN)要求编译器仅应用符合 IEEE-754 规则的优化。它可能不会将 x*0.0 转换为 0.0(这对于 -0、无穷大和 NaN 是不正确的)也可能不会优化 0.0-x 进入 -x，因为根据 IEEE-754 的第 5.5.1 节，否定是位级操作(对零和 NaN 产生不同的结果)。大多数编译器会提供一个标志，强制使用“安全”转换，例如-fp-model=precise 用于英特尔 C/C++ 编译器。

一个额外的警告适用于在参数缩减期间使用 nearbyint 函数。和rint一样，这个函数指定按照当前的舍入方式进行舍入。当未使用 fenv.h 时，舍入模式默认为“to-nearest-or-even”。使用它时，存在定向舍入模式生效的风险。这可以通过使用 round 来解决，它始终提供独立于当前舍入模式的舍入模式“舍入到最近，远离零”。然而，这个函数往往会变慢，因为它在大多数处理器架构上不受等效机器指令的支持。

关于性能的说明:下面的 C99 代码在很大程度上依赖于 fma() 的使用，它实现了 fused multiply-add。手术。在大多数现代硬件架构上，这由相应的硬件指令直接支持。如果不是这种情况，由于 FMA 仿真通常较慢，代码可能会显着变慢。

 #include <math.h>
 #include <stdint.h>

/* Writes result sine result sin(πa) to the location pointed to by sp
   Writes result cosine result cos(πa) to the location pointed to by cp

   In extensive testing, no errors > 0.97 ulp were found in either the sine
   or cosine results, suggesting the results returned are faithfully rounded.
*/
void my_sincospi (double a, double *sp, double *cp)
{
    double c, r, s, t, az;
    int64_t i;

    az = a * 0.0; // must be evaluated with IEEE-754 semantics
    /* for |a| >= 2**53, cospi(a) = 1.0, but cospi(Inf) = NaN */
    a = (fabs (a) < 9.0071992547409920e+15) ? a : az;  // 0x1.0p53
    /* reduce argument to primary approximation interval (-0.25, 0.25) */
    r = nearbyint (a + a); // must use IEEE-754 "to nearest" rounding
    i = (int64_t)r;
    t = fma (-0.5, r, a);
    /* compute core approximations */
    s = t * t;
    /* Approximate cos(pi*x) for x in [-0.25,0.25] */
    r =            -1.0369917389758117e-4;
    r = fma (r, s,  1.9294935641298806e-3);
    r = fma (r, s, -2.5806887942825395e-2);
    r = fma (r, s,  2.3533063028328211e-1);
    r = fma (r, s, -1.3352627688538006e+0);
    r = fma (r, s,  4.0587121264167623e+0);
    r = fma (r, s, -4.9348022005446790e+0);
    c = fma (r, s,  1.0000000000000000e+0);
    /* Approximate sin(pi*x) for x in [-0.25,0.25] */
    r =             4.6151442520157035e-4;
    r = fma (r, s, -7.3700183130883555e-3);
    r = fma (r, s,  8.2145868949323936e-2);
    r = fma (r, s, -5.9926452893214921e-1);
    r = fma (r, s,  2.5501640398732688e+0);
    r = fma (r, s, -5.1677127800499516e+0);
    s = s * t;
    r = r * s;
    s = fma (t, 3.1415926535897931e+0, r);
    /* map results according to quadrant */
    if (i & 2) {
        s = 0.0 - s; // must be evaluated with IEEE-754 semantics
        c = 0.0 - c; // must be evaluated with IEEE-754 semantics
    }
    if (i & 1) { 
        t = 0.0 - s; // must be evaluated with IEEE-754 semantics
        s = c;
        c = t;
    }
    /* IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n */
    if (a == floor (a)) s = az;
    *sp = s;
    *cp = c;
}

单精度版本基本上仅在核心近似方面有所不同。使用详尽测试可以精确确定误差范围。

#include <math.h>
#include <stdint.h>

/* Writes result sine result sin(πa) to the location pointed to by sp
   Writes result cosine result cos(πa) to the location pointed to by cp

   In exhaustive testing, the maximum error in sine results was 0.96677 ulp,
   the maximum error in cosine results was 0.96563 ulp, meaning results are
   faithfully rounded.
*/
void my_sincospif (float a, float *sp, float *cp)
{
    float az, t, c, r, s;
    int32_t i;

    az = a * 0.0f; // must be evaluated with IEEE-754 semantics
    /* for |a| > 2**24, cospi(a) = 1.0f, but cospi(Inf) = NaN */
    a = (fabsf (a) < 0x1.0p24f) ? a : az;
    r = nearbyintf (a + a); // must use IEEE-754 "to nearest" rounding
    i = (int32_t)r;
    t = fmaf (-0.5f, r, a);
    /* compute core approximations */
    s = t * t;
    /* Approximate cos(pi*x) for x in [-0.25,0.25] */
    r =              0x1.d9e000p-3f;
    r = fmaf (r, s, -0x1.55c400p+0f);
    r = fmaf (r, s,  0x1.03c1cep+2f);
    r = fmaf (r, s, -0x1.3bd3ccp+2f);
    c = fmaf (r, s,  0x1.000000p+0f);
    /* Approximate sin(pi*x) for x in [-0.25,0.25] */
    r =             -0x1.310000p-1f;
    r = fmaf (r, s,  0x1.46737ep+1f);
    r = fmaf (r, s, -0x1.4abbfep+2f);
    r = (t * s) * r;
    s = fmaf (t, 0x1.921fb6p+1f, r);
    if (i & 2) {
        s = 0.0f - s; // must be evaluated with IEEE-754 semantics
        c = 0.0f - c; // must be evaluated with IEEE-754 semantics
    }
    if (i & 1) {
        t = 0.0f - s; // must be evaluated with IEEE-754 semantics
        s = c;
        c = t;
    }
    /* IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n */
    if (a == floorf (a)) s = az;
    *sp = s;
    *cp = c;
}