diff options
Diffstat (limited to 'libm/ldouble/sinl.c')
| -rw-r--r-- | libm/ldouble/sinl.c | 342 |
1 files changed, 342 insertions, 0 deletions
diff --git a/libm/ldouble/sinl.c b/libm/ldouble/sinl.c new file mode 100644 index 000000000..dc7d739f9 --- /dev/null +++ b/libm/ldouble/sinl.c @@ -0,0 +1,342 @@ +/* sinl.c + * + * Circular sine, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, sinl(); + * + * y = sinl( x ); + * + * + * + * DESCRIPTION: + * + * Range reduction is into intervals of pi/4. The reduction + * error is nearly eliminated by contriving an extended precision + * modular arithmetic. + * + * Two polynomial approximating functions are employed. + * Between 0 and pi/4 the sine is approximated by the Cody + * and Waite polynomial form + * x + x**3 P(x**2) . + * Between pi/4 and pi/2 the cosine is represented as + * 1 - .5 x**2 + x**4 Q(x**2) . + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-5.5e11 200,000 1.2e-19 2.9e-20 + * + * ERROR MESSAGES: + * + * message condition value returned + * sin total loss x > 2**39 0.0 + * + * Loss of precision occurs for x > 2**39 = 5.49755813888e11. + * The routine as implemented flags a TLOSS error for + * x > 2**39 and returns 0.0. + */ + /* cosl.c + * + * Circular cosine, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, cosl(); + * + * y = cosl( x ); + * + * + * + * DESCRIPTION: + * + * Range reduction is into intervals of pi/4. The reduction + * error is nearly eliminated by contriving an extended precision + * modular arithmetic. + * + * Two polynomial approximating functions are employed. + * Between 0 and pi/4 the cosine is approximated by + * 1 - .5 x**2 + x**4 Q(x**2) . + * Between pi/4 and pi/2 the sine is represented by the Cody + * and Waite polynomial form + * x + x**3 P(x**2) . + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-5.5e11 50000 1.2e-19 2.9e-20 + */ + +/* sin.c */ + +/* +Cephes Math Library Release 2.7: May, 1998 +Copyright 1985, 1990, 1998 by Stephen L. Moshier +*/ + +#include <math.h> + +#ifdef UNK +static long double sincof[7] = { +-7.5785404094842805756289E-13L, + 1.6058363167320443249231E-10L, +-2.5052104881870868784055E-8L, + 2.7557319214064922217861E-6L, +-1.9841269841254799668344E-4L, + 8.3333333333333225058715E-3L, +-1.6666666666666666640255E-1L, +}; +static long double coscof[7] = { + 4.7377507964246204691685E-14L, +-1.1470284843425359765671E-11L, + 2.0876754287081521758361E-9L, +-2.7557319214999787979814E-7L, + 2.4801587301570552304991E-5L, +-1.3888888888888872993737E-3L, + 4.1666666666666666609054E-2L, +}; +static long double DP1 = 7.853981554508209228515625E-1L; +static long double DP2 = 7.946627356147928367136046290398E-9L; +static long double DP3 = 3.061616997868382943065164830688E-17L; +#endif + +#ifdef IBMPC +static short sincof[] = { +0x4e27,0xe1d6,0x2389,0xd551,0xbfd6, XPD +0x64d7,0xe706,0x4623,0xb090,0x3fde, XPD +0x01b1,0xbf34,0x2946,0xd732,0xbfe5, XPD +0xc8f7,0x9845,0x1d29,0xb8ef,0x3fec, XPD +0x6514,0x0c53,0x00d0,0xd00d,0xbff2, XPD +0x569a,0x8888,0x8888,0x8888,0x3ff8, XPD +0xaa97,0xaaaa,0xaaaa,0xaaaa,0xbffc, XPD +}; +static short coscof[] = { +0x7436,0x6f99,0x8c3a,0xd55e,0x3fd2, XPD +0x2f37,0x58f4,0x920f,0xc9c9,0xbfda, XPD +0x5350,0x659e,0xc648,0x8f76,0x3fe2, XPD +0x4d2b,0xf5c6,0x7dba,0x93f2,0xbfe9, XPD +0x53ed,0x0c66,0x00d0,0xd00d,0x3fef, XPD +0x7b67,0x0b60,0x60b6,0xb60b,0xbff5, XPD +0xaa9a,0xaaaa,0xaaaa,0xaaaa,0x3ffa, XPD +}; +static short P1[] = {0x0000,0x0000,0xda80,0xc90f,0x3ffe, XPD}; +static short P2[] = {0x0000,0x0000,0xa300,0x8885,0x3fe4, XPD}; +static short P3[] = {0x3707,0xa2e0,0x3198,0x8d31,0x3fc8, XPD}; +#define DP1 *(long double *)P1 +#define DP2 *(long double *)P2 +#define DP3 *(long double *)P3 +#endif + +#ifdef MIEEE +static long sincof[] = { +0xbfd60000,0xd5512389,0xe1d64e27, +0x3fde0000,0xb0904623,0xe70664d7, +0xbfe50000,0xd7322946,0xbf3401b1, +0x3fec0000,0xb8ef1d29,0x9845c8f7, +0xbff20000,0xd00d00d0,0x0c536514, +0x3ff80000,0x88888888,0x8888569a, +0xbffc0000,0xaaaaaaaa,0xaaaaaa97, +}; +static long coscof[] = { +0x3fd20000,0xd55e8c3a,0x6f997436, +0xbfda0000,0xc9c9920f,0x58f42f37, +0x3fe20000,0x8f76c648,0x659e5350, +0xbfe90000,0x93f27dba,0xf5c64d2b, +0x3fef0000,0xd00d00d0,0x0c6653ed, +0xbff50000,0xb60b60b6,0x0b607b67, +0x3ffa0000,0xaaaaaaaa,0xaaaaaa9a, +}; +static long P1[] = {0x3ffe0000,0xc90fda80,0x00000000}; +static long P2[] = {0x3fe40000,0x8885a300,0x00000000}; +static long P3[] = {0x3fc80000,0x8d313198,0xa2e03707}; +#define DP1 *(long double *)P1 +#define DP2 *(long double *)P2 +#define DP3 *(long double *)P3 +#endif + +static long double lossth = 5.49755813888e11L; /* 2^39 */ +extern long double PIO4L; +#ifdef ANSIPROT +extern long double polevll ( long double, void *, int ); +extern long double floorl ( long double ); +extern long double ldexpl ( long double, int ); +extern int isnanl ( long double ); +extern int isfinitel ( long double ); +#else +long double polevll(), floorl(), ldexpl(), isnanl(), isfinitel(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif + +long double sinl(x) +long double x; +{ +long double y, z, zz; +int j, sign; + +#ifdef NANS +if( isnanl(x) ) + return(x); +#endif +#ifdef MINUSZERO +if( x == 0.0L ) + return(x); +#endif +#ifdef NANS +if( !isfinitel(x) ) + { + mtherr( "sinl", DOMAIN ); +#ifdef NANS + return(NANL); +#else + return(0.0L); +#endif + } +#endif +/* make argument positive but save the sign */ +sign = 1; +if( x < 0 ) + { + x = -x; + sign = -1; + } + +if( x > lossth ) + { + mtherr( "sinl", TLOSS ); + return(0.0L); + } + +y = floorl( x/PIO4L ); /* integer part of x/PIO4 */ + +/* strip high bits of integer part to prevent integer overflow */ +z = ldexpl( y, -4 ); +z = floorl(z); /* integer part of y/8 */ +z = y - ldexpl( z, 4 ); /* y - 16 * (y/16) */ + +j = z; /* convert to integer for tests on the phase angle */ +/* map zeros to origin */ +if( j & 1 ) + { + j += 1; + y += 1.0L; + } +j = j & 07; /* octant modulo 360 degrees */ +/* reflect in x axis */ +if( j > 3) + { + sign = -sign; + j -= 4; + } + +/* Extended precision modular arithmetic */ +z = ((x - y * DP1) - y * DP2) - y * DP3; + +zz = z * z; +if( (j==1) || (j==2) ) + { + y = 1.0L - ldexpl(zz,-1) + zz * zz * polevll( zz, coscof, 6 ); + } +else + { + y = z + z * (zz * polevll( zz, sincof, 6 )); + } + +if(sign < 0) + y = -y; + +return(y); +} + + + + + +long double cosl(x) +long double x; +{ +long double y, z, zz; +long i; +int j, sign; + + +#ifdef NANS +if( isnanl(x) ) + return(x); +#endif +#ifdef INFINITIES +if( !isfinitel(x) ) + { + mtherr( "cosl", DOMAIN ); +#ifdef NANS + return(NANL); +#else + return(0.0L); +#endif + } +#endif + +/* make argument positive */ +sign = 1; +if( x < 0 ) + x = -x; + +if( x > lossth ) + { + mtherr( "cosl", TLOSS ); + return(0.0L); + } + +y = floorl( x/PIO4L ); +z = ldexpl( y, -4 ); +z = floorl(z); /* integer part of y/8 */ +z = y - ldexpl( z, 4 ); /* y - 16 * (y/16) */ + +/* integer and fractional part modulo one octant */ +i = z; +if( i & 1 ) /* map zeros to origin */ + { + i += 1; + y += 1.0L; + } +j = i & 07; +if( j > 3) + { + j -=4; + sign = -sign; + } + +if( j > 1 ) + sign = -sign; + +/* Extended precision modular arithmetic */ +z = ((x - y * DP1) - y * DP2) - y * DP3; + +zz = z * z; +if( (j==1) || (j==2) ) + { + y = z + z * (zz * polevll( zz, sincof, 6 )); + } +else + { + y = 1.0L - ldexpl(zz,-1) + zz * zz * polevll( zz, coscof, 6 ); + } + +if(sign < 0) + y = -y; + +return(y); +} |