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

Question

函数sinpi(x)计算sin(πx)，函数cospi(x)计算cos(πx)，其中与π的乘法隐含在函数中。这些函数最初是作为Sun Microsystems在20世纪80年代末中的扩展引入到C标准数学库中的。IEEE 754 sinPi -2008在第9节中指定了等效函数cosPi和™。

sin(πx)和cos(πx)是自然存在的。一个非常简单的例子是Box-Muller变换(G.E.P.Box和MervinE.Muller，“关于随机正态偏差产生的注记”)。“数理统计年鉴”，第29卷，第2期，第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功能作为系统库中的扩展。例如，苹果将__sinpi、__cospi和相应的单精度版本__sinpif、__cospif添加到iOS 7和OSX10.9 (介绍性，幻灯片101)中。但是对于许多其他平台来说，C程序没有现成的实现。

与传统的使用sin (M_PI * x)和cos (M_PI * x)的方法相比，sinpi和cospi的使用通过与π的内部乘法来减少舍入误差，从而提高了精度，同时由于参数的简化也带来了性能优势。

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

Answer 1

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

基本的算法结构很简单。由于非常大的参数总是偶数整数，因此2π的倍数，它们的正弦和余弦值是众所周知的。其他参数被折叠到范围-¼，+1/4，同时记录象限信息。利用多项式极小极大逼近计算初等逼近区间上的正弦和余弦。最后，通过周期性的结果交换和符号变化，利用象限数据将初步结果映射到最终结果。

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

还有一个附加的警告适用于在参数约简期间使用nearbyint函数。与rint一样，此函数根据当前的舍入模式指定为整数。当不使用fenv.h时，舍入模式默认为整“to -最近-或-偶数”。当它被使用时，有一个风险，一个定向舍入模式是有效的。这可以通过使用round来解决，它总是提供舍入模式“从零到最近，与当前舍入模式无关”。但是，由于在大多数处理器体系结构上不支持等效的机器指令，此函数将趋于较慢。

关于性能的注意事项:下面的C99代码严重依赖于fma()的使用，后者实现了融合相加操作。在大多数现代硬件体系结构中，这是由相应的硬件指令直接支持的。如果不是这样的话，代码可能会由于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;
}