diff options
Diffstat (limited to 'libm/ldouble/sinl.c')
| -rw-r--r-- | libm/ldouble/sinl.c | 342 |
1 files changed, 0 insertions, 342 deletions
diff --git a/libm/ldouble/sinl.c b/libm/ldouble/sinl.c deleted file mode 100644 index dc7d739f9..000000000 --- a/libm/ldouble/sinl.c +++ /dev/null @@ -1,342 +0,0 @@ -/* 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); -} |