/* sindg.c * * Circular sine of angle in degrees * * * * SYNOPSIS: * * double x, y, sindg(); * * y = sindg( x ); * * * * DESCRIPTION: * * Range reduction is into intervals of 45 degrees. * * Two polynomial approximating functions are employed. * Between 0 and pi/4 the sine is approximated by * x + x**3 P(x**2). * Between pi/4 and pi/2 the cosine is represented as * 1 - x**2 P(x**2). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * DEC +-1000 3100 3.3e-17 9.0e-18 * IEEE +-1000 30000 2.3e-16 5.6e-17 * * ERROR MESSAGES: * * message condition value returned * sindg total loss x > 8.0e14 (DEC) 0.0 * x > 1.0e14 (IEEE) * */ /* cosdg.c * * Circular cosine of angle in degrees * * * * SYNOPSIS: * * double x, y, cosdg(); * * y = cosdg( x ); * * * * DESCRIPTION: * * Range reduction is into intervals of 45 degrees. * * Two polynomial approximating functions are employed. * Between 0 and pi/4 the cosine is approximated by * 1 - x**2 P(x**2). * Between pi/4 and pi/2 the sine is represented as * x + x**3 P(x**2). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * DEC +-1000 3400 3.5e-17 9.1e-18 * IEEE +-1000 30000 2.1e-16 5.7e-17 * See also sin(). * */ /* Cephes Math Library Release 2.0: April, 1987 * Copyright 1985, 1987 by Stephen L. Moshier * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ #include "mconf.h" #ifdef UNK static double sincof[] = { 1.58962301572218447952E-10, -2.50507477628503540135E-8, 2.75573136213856773549E-6, -1.98412698295895384658E-4, 8.33333333332211858862E-3, -1.66666666666666307295E-1}; static double coscof[] = { 1.13678171382044553091E-11, -2.08758833757683644217E-9, 2.75573155429816611547E-7, -2.48015872936186303776E-5, 1.38888888888806666760E-3, -4.16666666666666348141E-2, 4.99999999999999999798E-1}; static double PI180 = 1.74532925199432957692E-2; /* pi/180 */ static double lossth = 1.0e14; #endif #ifdef DEC static unsigned short sincof[] = { 0030056, 0143750, 0177170, 0073013, 0131727, 0027455, 0044510, 0132205, 0033470, 0167432, 0131752, 0042263, 0135120, 0006400, 0146776, 0174027, 0036410, 0104210, 0104207, 0137202, 0137452, 0125252, 0125252, 0125103}; static unsigned short coscof[] = { 0027107, 0176030, 0153315, 0110312, 0131017, 0072476, 0007450, 0123243, 0032623, 0171174, 0070066, 0146445, 0134320, 0006400, 0147355, 0163313, 0035666, 0005540, 0133012, 0165067, 0137052, 0125252, 0125252, 0125206, 0040000, 0000000, 0000000, 0000000}; static unsigned short P1[] = {0036616, 0175065, 0011224, 0164711}; #define PI180 *(double *)P1 static double lossth = 8.0e14; #endif #ifdef IBMPC static unsigned short sincof[] = { 0x0ec1, 0x1fcf, 0xd8fd, 0x3de5, 0x1691, 0xa929, 0xe5e5, 0xbe5a, 0x4896, 0x567d, 0x1de3, 0x3ec7, 0xdf03, 0x19bf, 0x01a0, 0xbf2a, 0xf7d0, 0x1110, 0x1111, 0x3f81, 0x5548, 0x5555, 0x5555, 0xbfc5}; static unsigned short coscof[] = { 0xb219, 0x1ad9, 0xff83, 0x3da8, 0x14d4, 0xc1e5, 0xeea7, 0xbe21, 0xd9a5, 0x8e06, 0x7e4f, 0x3e92, 0xbcd9, 0x19dd, 0x01a0, 0xbefa, 0x5d47, 0x16c1, 0xc16c, 0x3f56, 0x5551, 0x5555, 0x5555, 0xbfa5, 0x0000, 0x0000, 0x0000, 0x3fe0}; static unsigned short P1[] = {0x9d39, 0xa252, 0xdf46, 0x3f91}; #define PI180 *(double *)P1 static double lossth = 1.0e14; #endif #ifdef MIEEE static unsigned short sincof[] = { 0x3de5, 0xd8fd, 0x1fcf, 0x0ec1, 0xbe5a, 0xe5e5, 0xa929, 0x1691, 0x3ec7, 0x1de3, 0x567d, 0x4896, 0xbf2a, 0x01a0, 0x19bf, 0xdf03, 0x3f81, 0x1111, 0x1110, 0xf7d0, 0xbfc5, 0x5555, 0x5555, 0x5548}; static unsigned short coscof[] = { 0x3da8, 0xff83, 0x1ad9, 0xb219, 0xbe21, 0xeea7, 0xc1e5, 0x14d4, 0x3e92, 0x7e4f, 0x8e06, 0xd9a5, 0xbefa, 0x01a0, 0x19dd, 0xbcd9, 0x3f56, 0xc16c, 0x16c1, 0x5d47, 0xbfa5, 0x5555, 0x5555, 0x5551, 0x3fe0, 0x0000, 0x0000, 0x0000}; static unsigned short P1[] = {0x3f91, 0xdf46, 0xa252, 0x9d39}; #define PI180 *(double *)P1 static double lossth = 1.0e14; #endif #ifdef ANSIPROT extern double polevl(double, void *, int); extern double floor(double); extern double ldexp(double, int); #else double polevl(), floor(), ldexp(); #endif extern double PIO4; double sindg(x) double x; { double y, z, zz; int j, sign; /* make argument positive but save the sign */ sign = 1; if (x < 0) { x = -x; sign = -1; } if (x > lossth) { mtherr("sindg", TLOSS); return (0.0); } y = floor(x / 45.0); /* integer part of x/PIO4 */ /* strip high bits of integer part to prevent integer overflow */ z = ldexp(y, -4); z = floor(z); /* integer part of y/8 */ z = y - ldexp(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.0; } j = j & 07; /* octant modulo 360 degrees */ /* reflect in x axis */ if (j > 3) { sign = -sign; j -= 4; } z = x - y * 45.0; /* x mod 45 degrees */ z *= PI180; /* multiply by pi/180 to convert to radians */ zz = z * z; if ((j == 1) || (j == 2)) { y = 1.0 - zz * polevl(zz, coscof, 6); } else { y = z + z * (zz * polevl(zz, sincof, 5)); } if (sign < 0) y = -y; return (y); } double cosdg(x) double x; { double y, z, zz; int j, sign; /* make argument positive */ sign = 1; if (x < 0) x = -x; if (x > lossth) { mtherr("cosdg", TLOSS); return (0.0); } y = floor(x / 45.0); z = ldexp(y, -4); z = floor(z); /* integer part of y/8 */ z = y - ldexp(z, 4); /* y - 16 * (y/16) */ /* integer and fractional part modulo one octant */ j = z; if (j & 1) /* map zeros to origin */ { j += 1; y += 1.0; } j = j & 07; if (j > 3) { j -= 4; sign = -sign; } if (j > 1) sign = -sign; z = x - y * 45.0; /* x mod 45 degrees */ z *= PI180; /* multiply by pi/180 to convert to radians */ zz = z * z; if ((j == 1) || (j == 2)) { y = z + z * (zz * polevl(zz, sincof, 5)); } else { y = 1.0 - zz * polevl(zz, coscof, 6); } if (sign < 0) y = -y; return (y); }