https://gcc.gnu.org/bugzilla/show_bug.cgi?id=116299
--- Comment #2 from Ryan <rmaguire314 at gmail dot com> ---
It is a part of a rather large library.
How much of it do you want?
Here is the full function, in a working condition.
/* TMPL_STATIC_INLINE macro found here.
*/
#include <libtmpl/include/tmpl_config.h>
/* The denominator of the asymptotic expansion is scaled by pi.
*/
#define TMPL_ONE_PI (+3.14159265358979323846264338327950288419716939E+00L)
/* Used to compute sin(pi t) and cos(pi t) simultaneously.
*/
extern void
tmpl_LDouble_SinCosPi(long double t, long double *sin_t, long double *cos_t);
/* Different splitting values are needed, depending on the type. If the *
* mantissa is N bits, the magic number is 2^(N - floor(N/3) + 1) + 1. This *
* guarantees that for |x| > 2^floor(N/3) the high part is even, meaning *
* xhi^2 / 2 is an integer. Since sin(pi x) and cos(pi x) are periodic, we *
* can then discard this computation entirely. We need only concentrate on *
* the lower N - floor(N/3) bits.
*/
#if TMPL_LDOUBLE_TYPE == TMPL_LDOUBLE_64_BIT
/* 52-bit mantissa, value is 2^36 + 1.
*/
#define TMPL_LDOUBLE_SPLIT (6.8719476737E+10L)
/* Quadruple precision has a much large mantissa, hence a larger value.
*/
#elif TMPL_LDOUBLE_TYPE == TMPL_LDOUBLE_128_BIT
/* 112-bit mantissa, value is 2^76 + 1.
*/
#define TMPL_LDOUBLE_SPLIT (+7.5557863725914323419137E+22L)
/* Double-double, special case, we split the higher double.
*/
#elif TMPL_LDOUBLE_TYPE == TMPL_LDOUBLE_DOUBLEDOUBLE
/* Double-double, split by 2^35 + 1 (splitting as a "double").
*/
#define TMPL_LDOUBLE_SPLIT (+3.4359738369E+10)
/* The most common version, 80-bit extended / portable.
*/
#else
/* 63-bit mantissa, value is 2^43 + 1.
*/
#define TMPL_LDOUBLE_SPLIT (+8.796093022209E+12L)
#endif
/* End of double vs. extended vs. double-double vs. quadruple.
*/
/* Function for computing the normalized Fresnel cosine of a large input.
*/
TMPL_STATIC_INLINE
long double tmpl_LDouble_Normalized_Fresnel_Cos_Asymptotic(long double x)
{
/* Use the double-double trick, split x into two parts, high and low.
*/
#if TMPL_LDOUBLE_TYPE == TMPL_LDOUBLE_DOUBLEDOUBLE
volatile const double x_double = (double)x;
volatile const double split = x_double*3.4359738369E+10;
volatile const double xhi_double = split - (split - x_double);
const long double xhi = (long double)xhi_double;
#else
const long double split = TMPL_LDOUBLE_SPLIT * x;
const long double xhi = split - (split - x);
#endif
const long double xlo = x - xhi;
/* The scale factor for the asymptotic expansion. For large x we only *
* need the first term of the approximation.
*/
const long double t = 1.0L / (TMPL_ONE_PI * x);
/* For large x we have xhi^2 / 2 is an even integer. Since sin(pi t) *
* is periodic with period 2, the xhi^2 term can be disregarded. The *
* argument we then care about is pi (2 xhi xlo + xlo^2) / 2.
*/
long double sin_hi, cos_hi, sin_lo, cos_lo, sin_x;
/* Compute sin(pi/2 (2 xhi xlo + xlo^2)) using the angle sum formula.
*/
tmpl_LDouble_SinCosPi(xlo * xhi, &sin_hi, &cos_hi);
tmpl_LDouble_SinCosPi(0.5L * xlo * xlo, &sin_lo, &cos_lo);
sin_x = cos_hi*sin_lo + cos_lo*sin_hi;
/* The first term of the asymptotic expansion is all that is needed.
*/
return 0.5L + t * sin_x;
}
