tan.c 6.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
  1. /* tan.c
  2. *
  3. * Circular tangent
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, tan();
  10. *
  11. * y = tan( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns the circular tangent of the radian argument x.
  18. *
  19. * Range reduction is modulo pi/4. A rational function
  20. * x + x**3 P(x**2)/Q(x**2)
  21. * is employed in the basic interval [0, pi/4].
  22. *
  23. *
  24. *
  25. * ACCURACY:
  26. *
  27. * Relative error:
  28. * arithmetic domain # trials peak rms
  29. * DEC +-1.07e9 44000 4.1e-17 1.0e-17
  30. * IEEE +-1.07e9 30000 2.9e-16 8.1e-17
  31. *
  32. * ERROR MESSAGES:
  33. *
  34. * message condition value returned
  35. * tan total loss x > 1.073741824e9 0.0
  36. *
  37. */
  38. /* cot.c
  39. *
  40. * Circular cotangent
  41. *
  42. *
  43. *
  44. * SYNOPSIS:
  45. *
  46. * double x, y, cot();
  47. *
  48. * y = cot( x );
  49. *
  50. *
  51. *
  52. * DESCRIPTION:
  53. *
  54. * Returns the circular cotangent of the radian argument x.
  55. *
  56. * Range reduction is modulo pi/4. A rational function
  57. * x + x**3 P(x**2)/Q(x**2)
  58. * is employed in the basic interval [0, pi/4].
  59. *
  60. *
  61. *
  62. * ACCURACY:
  63. *
  64. * Relative error:
  65. * arithmetic domain # trials peak rms
  66. * IEEE +-1.07e9 30000 2.9e-16 8.2e-17
  67. *
  68. *
  69. * ERROR MESSAGES:
  70. *
  71. * message condition value returned
  72. * cot total loss x > 1.073741824e9 0.0
  73. * cot singularity x = 0 INFINITY
  74. *
  75. */
  76. /*
  77. Cephes Math Library Release 2.8: June, 2000
  78. yright 1984, 1995, 2000 by Stephen L. Moshier
  79. */
  80. #include "mconf.h"
  81. #ifdef UNK
  82. static double P[] = {-1.30936939181383777646E4, 1.15351664838587416140E6,
  83. -1.79565251976484877988E7};
  84. static double Q[] = {
  85. /* 1.00000000000000000000E0,*/
  86. 1.36812963470692954678E4, -1.32089234440210967447E6,
  87. 2.50083801823357915839E7, -5.38695755929454629881E7};
  88. static double DP1 = 7.853981554508209228515625E-1;
  89. static double DP2 = 7.94662735614792836714E-9;
  90. static double DP3 = 3.06161699786838294307E-17;
  91. static double lossth = 1.073741824e9;
  92. #endif
  93. #ifdef DEC
  94. static unsigned short P[] = {0143514, 0113306, 0111171, 0174674,
  95. 0045214, 0147545, 0027744, 0167346,
  96. 0146210, 0177526, 0114514, 0105660};
  97. static unsigned short Q[] = {
  98. /*0040200,0000000,0000000,0000000,*/
  99. 0043525, 0142457, 0072633, 0025617, 0145241, 0036742, 0140525, 0162256,
  100. 0046276, 0146176, 0013526, 0143573, 0146515, 0077401, 0162762, 0150607};
  101. /* 7.853981629014015197753906250000E-1 */
  102. static unsigned short P1[] = {
  103. 0040111,
  104. 0007732,
  105. 0120000,
  106. 0000000,
  107. };
  108. /* 4.960467869796758577649598009884E-10 */
  109. static unsigned short P2[] = {
  110. 0030410,
  111. 0055060,
  112. 0100000,
  113. 0000000,
  114. };
  115. /* 2.860594363054915898381331279295E-18 */
  116. static unsigned short P3[] = {
  117. 0021523,
  118. 0011431,
  119. 0105056,
  120. 0001560,
  121. };
  122. #define DP1 *(double *)P1
  123. #define DP2 *(double *)P2
  124. #define DP3 *(double *)P3
  125. static double lossth = 1.073741824e9;
  126. #endif
  127. #ifdef IBMPC
  128. static unsigned short P[] = {0x3f38, 0xd24f, 0x92d8, 0xc0c9, 0x9ddd, 0xa5fc,
  129. 0x99ec, 0x4131, 0x9176, 0xd329, 0x1fea, 0xc171};
  130. static unsigned short Q[] = {
  131. /*0x0000,0x0000,0x0000,0x3ff0,*/
  132. 0x6572, 0xeeb3, 0xb8a5, 0x40ca, 0xbc96, 0x582a, 0x27bc, 0xc134,
  133. 0xd8ef, 0xc2ea, 0xd98f, 0x4177, 0x5a31, 0x3cbe, 0xafe0, 0xc189};
  134. /*
  135. 7.85398125648498535156E-1,
  136. 3.77489470793079817668E-8,
  137. 2.69515142907905952645E-15,
  138. */
  139. static unsigned short P1[] = {0x0000, 0x4000, 0x21fb, 0x3fe9};
  140. static unsigned short P2[] = {0x0000, 0x0000, 0x442d, 0x3e64};
  141. static unsigned short P3[] = {0x5170, 0x98cc, 0x4698, 0x3ce8};
  142. #define DP1 *(double *)P1
  143. #define DP2 *(double *)P2
  144. #define DP3 *(double *)P3
  145. static double lossth = 1.073741824e9;
  146. #endif
  147. #ifdef MIEEE
  148. static unsigned short P[] = {0xc0c9, 0x92d8, 0xd24f, 0x3f38, 0x4131, 0x99ec,
  149. 0xa5fc, 0x9ddd, 0xc171, 0x1fea, 0xd329, 0x9176};
  150. static unsigned short Q[] = {0x40ca, 0xb8a5, 0xeeb3, 0x6572, 0xc134, 0x27bc,
  151. 0x582a, 0xbc96, 0x4177, 0xd98f, 0xc2ea, 0xd8ef,
  152. 0xc189, 0xafe0, 0x3cbe, 0x5a31};
  153. static unsigned short P1[] = {0x3fe9, 0x21fb, 0x4000, 0x0000};
  154. static unsigned short P2[] = {0x3e64, 0x442d, 0x0000, 0x0000};
  155. static unsigned short P3[] = {
  156. 0x3ce8,
  157. 0x4698,
  158. 0x98cc,
  159. 0x5170,
  160. };
  161. #define DP1 *(double *)P1
  162. #define DP2 *(double *)P2
  163. #define DP3 *(double *)P3
  164. static double lossth = 1.073741824e9;
  165. #endif
  166. #ifdef ANSIPROT
  167. extern double polevl(double, void *, int);
  168. extern double p1evl(double, void *, int);
  169. extern double floor(double);
  170. extern double ldexp(double, int);
  171. extern int isnan(double);
  172. extern int isfinite(double);
  173. static double tancot(double, int);
  174. #else
  175. double polevl(), p1evl(), floor(), ldexp();
  176. static double tancot();
  177. int isnan(), isfinite();
  178. #endif
  179. extern double PIO4;
  180. extern double INFINITY;
  181. extern double NAN;
  182. double tan(x) double x;
  183. {
  184. #ifdef MINUSZERO
  185. if (x == 0.0)
  186. return (x);
  187. #endif
  188. #ifdef NANS
  189. if (isnan(x))
  190. return (x);
  191. if (!isfinite(x)) {
  192. mtherr("tan", DOMAIN);
  193. return (NAN);
  194. }
  195. #endif
  196. return (tancot(x, 0));
  197. }
  198. double cot(x) double x;
  199. {
  200. if (x == 0.0) {
  201. mtherr("cot", SING);
  202. return (INFINITY);
  203. }
  204. return (tancot(x, 1));
  205. }
  206. static double tancot(xx, cotflg) double xx;
  207. int cotflg;
  208. {
  209. double x, y, z, zz;
  210. int j, sign;
  211. /* make argument positive but save the sign */
  212. if (xx < 0) {
  213. x = -xx;
  214. sign = -1;
  215. } else {
  216. x = xx;
  217. sign = 1;
  218. }
  219. if (x > lossth) {
  220. if (cotflg)
  221. mtherr("cot", TLOSS);
  222. else
  223. mtherr("tan", TLOSS);
  224. return (0.0);
  225. }
  226. /* compute x mod PIO4 */
  227. y = floor(x / PIO4);
  228. /* strip high bits of integer part */
  229. z = ldexp(y, -3);
  230. z = floor(z); /* integer part of y/8 */
  231. z = y - ldexp(z, 3); /* y - 16 * (y/16) */
  232. /* integer and fractional part modulo one octant */
  233. j = z;
  234. /* map zeros and singularities to origin */
  235. if (j & 1) {
  236. j += 1;
  237. y += 1.0;
  238. }
  239. z = ((x - y * DP1) - y * DP2) - y * DP3;
  240. zz = z * z;
  241. if (zz > 1.0e-14)
  242. y = z + z * (zz * polevl(zz, P, 2) / p1evl(zz, Q, 4));
  243. else
  244. y = z;
  245. if (j & 2) {
  246. if (cotflg)
  247. y = -y;
  248. else
  249. y = -1.0 / y;
  250. } else {
  251. if (cotflg)
  252. y = 1.0 / y;
  253. }
  254. if (sign < 0)
  255. y = -y;
  256. return (y);
  257. }