使用标准C数学库实现sinpi（）和cospi（）
函数sinpi(x)计算sin（πx），函数cospi(x)计算cos（πx），其中与π的乘法隐含在函数内部。 这些函数最初被引入C标准数学库，作为Sun Microsystems在20世纪80年代后期的扩展。 IEEE Std 754™-2008在第9节中规定了等效函数sinPi和cosPi 。
有许多计算，sin（πx）和cos（πx）自然发生。 一个非常简单的例子是Box-Muller变换（GEP Box和Mervin E. Muller，“关于随机正态偏差的一个注记”。 “数学统计年鉴” ，第29卷，第2期，第610-611页） ），给定两个独立的随机变量U 1和U 2均匀分布，产生具有标准正态分布的独立随机变量Z 1和Z 2：
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, eg 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和cospifunction作为系统库中的扩展。 例如，Apple将__sinpi ， __cospi和相应的单精度版本__sinpif ， __cospif到iOS 7和OS X 10.9（ 演示文稿 ，幻灯片101）。 但对于许多其他平台，C程序无法轻松访问。
与使用例如sin (M_PI * x)和cos (M_PI * x)的传统方法相比， sinpi和cospi的使用通过使用π的内部乘法减少舍入误差来提高精度，并且还提供了性能优势更简单的参数减少。
如何使用标准C数学库以合理有效和标准兼容的方式实现sinpi()和cospi()function？
为简单起见，我将重点介绍sincospi() ，它同时提供正弦和余弦结果。 然后可以将sinpi和cospi构造为丢弃不需要的数据的包装函数。 在许多应用程序中，不需要处理浮点标志（请参阅fenv.h ），也不需要在大多数情况下报告错误报告，因此我将省略这些。
基本的算法结构很简单。 因为非常大的参数总是偶数整数，因此因此是2π的倍数，它们的正弦和余弦值是众所周知的。 在记录象限信息时，其他参数被折叠到[-¼，+¼]范围内。 多项式极小极大近似用于计算初级近似区间的正弦和余弦。 最后，象限数据用于通过周期性交换结果和符号变化将初步结果映射到最终结果。
正确处理特殊操作数（特别是-0，无穷大和NaN）需要编译器仅应用符合IEEE-754规则的优化。 它可能不会将x*0.0转换为0.0 （这对于-0，无穷大和NaN不正确）也不能将0.0-x优化为-x因为根据IEEE的5.5.1节，否定是一个位级操作754（对零和NaN产生不同的结果）。 大多数编译器都会提供一个强制使用“安全”转换的标志，例如英特尔C / C ++编译器的-fp-model=precise 。
另一个警告适用于在参数减少期间使用nearbyint函数。 与rint一样，此函数指定为根据当前舍入模式进行舍入。 如果不使用fenv.h ，则舍入模式默认为“to-nearest-or-even”。 使用时，存在定向舍入模式有效的风险。 这可以通过使用round ，它总是提供舍入模式“舍入到最近，从零开始”，与当前舍入模式无关。 但是，此function往往较慢，因为大多数处理器体系结构上的等效机器指令都不支持此function。
关于性能的说明：下面的C99代码在很大程度上依赖于fma()的使用， fma()实现了融合的乘法 – 加法运算。 在大多数现代硬件架构中，这由相应的硬件指令直接支持。 在不是这种情况的情况下，由于通常较慢的FMA仿真，代码可能经历显着的减速。
#include #include /* 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 #include /* 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; }