| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280 |
- /* tan.c
- *
- * Circular tangent
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, tan();
- *
- * y = tan( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns the circular tangent of the radian argument x.
- *
- * Range reduction is modulo pi/4. A rational function
- * x + x**3 P(x**2)/Q(x**2)
- * is employed in the basic interval [0, pi/4].
- *
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * DEC +-1.07e9 44000 4.1e-17 1.0e-17
- * IEEE +-1.07e9 30000 2.9e-16 8.1e-17
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * tan total loss x > 1.073741824e9 0.0
- *
- */
- /* cot.c
- *
- * Circular cotangent
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, cot();
- *
- * y = cot( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns the circular cotangent of the radian argument x.
- *
- * Range reduction is modulo pi/4. A rational function
- * x + x**3 P(x**2)/Q(x**2)
- * is employed in the basic interval [0, pi/4].
- *
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE +-1.07e9 30000 2.9e-16 8.2e-17
- *
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * cot total loss x > 1.073741824e9 0.0
- * cot singularity x = 0 INFINITY
- *
- */
- /*
- Cephes Math Library Release 2.8: June, 2000
- yright 1984, 1995, 2000 by Stephen L. Moshier
- */
- #include "mconf.h"
- #ifdef UNK
- static double P[] = {-1.30936939181383777646E4, 1.15351664838587416140E6,
- -1.79565251976484877988E7};
- static double Q[] = {
- /* 1.00000000000000000000E0,*/
- 1.36812963470692954678E4, -1.32089234440210967447E6,
- 2.50083801823357915839E7, -5.38695755929454629881E7};
- static double DP1 = 7.853981554508209228515625E-1;
- static double DP2 = 7.94662735614792836714E-9;
- static double DP3 = 3.06161699786838294307E-17;
- static double lossth = 1.073741824e9;
- #endif
- #ifdef DEC
- static unsigned short P[] = {0143514, 0113306, 0111171, 0174674,
- 0045214, 0147545, 0027744, 0167346,
- 0146210, 0177526, 0114514, 0105660};
- static unsigned short Q[] = {
- /*0040200,0000000,0000000,0000000,*/
- 0043525, 0142457, 0072633, 0025617, 0145241, 0036742, 0140525, 0162256,
- 0046276, 0146176, 0013526, 0143573, 0146515, 0077401, 0162762, 0150607};
- /* 7.853981629014015197753906250000E-1 */
- static unsigned short P1[] = {
- 0040111,
- 0007732,
- 0120000,
- 0000000,
- };
- /* 4.960467869796758577649598009884E-10 */
- static unsigned short P2[] = {
- 0030410,
- 0055060,
- 0100000,
- 0000000,
- };
- /* 2.860594363054915898381331279295E-18 */
- static unsigned short P3[] = {
- 0021523,
- 0011431,
- 0105056,
- 0001560,
- };
- #define DP1 *(double *)P1
- #define DP2 *(double *)P2
- #define DP3 *(double *)P3
- static double lossth = 1.073741824e9;
- #endif
- #ifdef IBMPC
- static unsigned short P[] = {0x3f38, 0xd24f, 0x92d8, 0xc0c9, 0x9ddd, 0xa5fc,
- 0x99ec, 0x4131, 0x9176, 0xd329, 0x1fea, 0xc171};
- static unsigned short Q[] = {
- /*0x0000,0x0000,0x0000,0x3ff0,*/
- 0x6572, 0xeeb3, 0xb8a5, 0x40ca, 0xbc96, 0x582a, 0x27bc, 0xc134,
- 0xd8ef, 0xc2ea, 0xd98f, 0x4177, 0x5a31, 0x3cbe, 0xafe0, 0xc189};
- /*
- 7.85398125648498535156E-1,
- 3.77489470793079817668E-8,
- 2.69515142907905952645E-15,
- */
- static unsigned short P1[] = {0x0000, 0x4000, 0x21fb, 0x3fe9};
- static unsigned short P2[] = {0x0000, 0x0000, 0x442d, 0x3e64};
- static unsigned short P3[] = {0x5170, 0x98cc, 0x4698, 0x3ce8};
- #define DP1 *(double *)P1
- #define DP2 *(double *)P2
- #define DP3 *(double *)P3
- static double lossth = 1.073741824e9;
- #endif
- #ifdef MIEEE
- static unsigned short P[] = {0xc0c9, 0x92d8, 0xd24f, 0x3f38, 0x4131, 0x99ec,
- 0xa5fc, 0x9ddd, 0xc171, 0x1fea, 0xd329, 0x9176};
- static unsigned short Q[] = {0x40ca, 0xb8a5, 0xeeb3, 0x6572, 0xc134, 0x27bc,
- 0x582a, 0xbc96, 0x4177, 0xd98f, 0xc2ea, 0xd8ef,
- 0xc189, 0xafe0, 0x3cbe, 0x5a31};
- static unsigned short P1[] = {0x3fe9, 0x21fb, 0x4000, 0x0000};
- static unsigned short P2[] = {0x3e64, 0x442d, 0x0000, 0x0000};
- static unsigned short P3[] = {
- 0x3ce8,
- 0x4698,
- 0x98cc,
- 0x5170,
- };
- #define DP1 *(double *)P1
- #define DP2 *(double *)P2
- #define DP3 *(double *)P3
- static double lossth = 1.073741824e9;
- #endif
- #ifdef ANSIPROT
- extern double polevl(double, void *, int);
- extern double p1evl(double, void *, int);
- extern double floor(double);
- extern double ldexp(double, int);
- extern int isnan(double);
- extern int isfinite(double);
- static double tancot(double, int);
- #else
- double polevl(), p1evl(), floor(), ldexp();
- static double tancot();
- int isnan(), isfinite();
- #endif
- extern double PIO4;
- extern double INFINITY;
- extern double NAN;
- double tan(x) double x;
- {
- #ifdef MINUSZERO
- if (x == 0.0)
- return (x);
- #endif
- #ifdef NANS
- if (isnan(x))
- return (x);
- if (!isfinite(x)) {
- mtherr("tan", DOMAIN);
- return (NAN);
- }
- #endif
- return (tancot(x, 0));
- }
- double cot(x) double x;
- {
- if (x == 0.0) {
- mtherr("cot", SING);
- return (INFINITY);
- }
- return (tancot(x, 1));
- }
- static double tancot(xx, cotflg) double xx;
- int cotflg;
- {
- double x, y, z, zz;
- int j, sign;
- /* make argument positive but save the sign */
- if (xx < 0) {
- x = -xx;
- sign = -1;
- } else {
- x = xx;
- sign = 1;
- }
- if (x > lossth) {
- if (cotflg)
- mtherr("cot", TLOSS);
- else
- mtherr("tan", TLOSS);
- return (0.0);
- }
- /* compute x mod PIO4 */
- y = floor(x / PIO4);
- /* strip high bits of integer part */
- z = ldexp(y, -3);
- z = floor(z); /* integer part of y/8 */
- z = y - ldexp(z, 3); /* y - 16 * (y/16) */
- /* integer and fractional part modulo one octant */
- j = z;
- /* map zeros and singularities to origin */
- if (j & 1) {
- j += 1;
- y += 1.0;
- }
- z = ((x - y * DP1) - y * DP2) - y * DP3;
- zz = z * z;
- if (zz > 1.0e-14)
- y = z + z * (zz * polevl(zz, P, 2) / p1evl(zz, Q, 4));
- else
- y = z;
- if (j & 2) {
- if (cotflg)
- y = -y;
- else
- y = -1.0 / y;
- } else {
- if (cotflg)
- y = 1.0 / y;
- }
- if (sign < 0)
- y = -y;
- return (y);
- }
|