The function sinpi(x)
computes sin(πx), and the function cospi(x)
computes cos(πx), where the multiplication with π is implicit inside the functions. These functions were initially introduced into the C standard math library as an extension by Sun Microsystems in the late 1980s. IEEE Std 754™-2008 specifies the equivalent functions sinPi
and cosPi
in section 9.
There are numerous computations where sin(πx) and cos(πx) occur naturally. A very simple example is the Box-Muller transform (G. E. P. Box and Mervin E. Muller, "A Note on the Generation of Random Normal Deviates". The Annals of Mathematical Statistics, Vol. 29, No. 2, pp. 610 - 611), which, given two independent random variables U₁ and U₂ with uniform distribution, produces independent random variables Z₁ and Z₂ with standard normal distribution:
Z₁ = √(-2 ln U₁) cos (2 π U₂)
Z₂ = √(-2 ln U₁) sin (2 π U₂)
A further example is the computation of sine and cosine for degree arguments, as in this computation of great-circle distance using the Haversine formula:
/* 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;
}
For C++, the Boost library provides sin_pi
and
cos_pi
, and some vendors offer sinpi
and cospi
functionality as extensions in system libraries. For example, Apple added __sinpi
, __cospi
and the corresponding single-precision versions __sinpif
, __cospif
to iOS 7 and OS X 10.9 (presentation, slide 101). But for many other platforms, there is no implementation readily accessible to C programs.
Compared with a traditional approach that uses e.g. sin (M_PI * x)
and cos (M_PI * x)
, the use of sinpi
and cospi
improves accuracy by reducing rounding error via the internal multiplication with π, and also offers performance advantages due to the much simpler argument reduction.
How can one use the standard C math library to implement sinpi()
and cospi()
functionality in a reasonably efficient and standard compliant fashion?
fenv()
orfesetround()
) to truncate/round-towards-zero is necessary. That way we can use e.g. Kahan sum/compensated sum, and split high-precision coefficients to several different limited-precision factors. Every other approach seems to rely on specific hardware (likefma()
, for which emulation is horribly slow) or implementation details. – Nominal Animal