123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256 |
- /* 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);
- }
|