tandg.c 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216
  1. /* tandg.c
  2. *
  3. * Circular tangent of argument in degrees
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, tandg();
  10. *
  11. * y = tandg( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns the circular tangent of the argument x in degrees.
  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 0,10 8000 3.4e-17 1.2e-17
  30. * IEEE 0,10 30000 3.2e-16 8.4e-17
  31. *
  32. * ERROR MESSAGES:
  33. *
  34. * message condition value returned
  35. * tandg total loss x > 8.0e14 (DEC) 0.0
  36. * x > 1.0e14 (IEEE)
  37. * tandg singularity x = 180 k + 90 MAXNUM
  38. */
  39. /* cotdg.c
  40. *
  41. * Circular cotangent of argument in degrees
  42. *
  43. *
  44. *
  45. * SYNOPSIS:
  46. *
  47. * double x, y, cotdg();
  48. *
  49. * y = cotdg( x );
  50. *
  51. *
  52. *
  53. * DESCRIPTION:
  54. *
  55. * Returns the circular cotangent of the argument x in degrees.
  56. *
  57. * Range reduction is modulo pi/4. A rational function
  58. * x + x**3 P(x**2)/Q(x**2)
  59. * is employed in the basic interval [0, pi/4].
  60. *
  61. *
  62. * ERROR MESSAGES:
  63. *
  64. * message condition value returned
  65. * cotdg total loss x > 8.0e14 (DEC) 0.0
  66. * x > 1.0e14 (IEEE)
  67. * cotdg singularity x = 180 k MAXNUM
  68. */
  69. /*
  70. Cephes Math Library Release 2.8: June, 2000
  71. Copyright 1984, 1987, 2000 by Stephen L. Moshier
  72. */
  73. #include "mconf.h"
  74. #ifdef UNK
  75. static double P[] = {-1.30936939181383777646E4, 1.15351664838587416140E6,
  76. -1.79565251976484877988E7};
  77. static double Q[] = {
  78. /* 1.00000000000000000000E0,*/
  79. 1.36812963470692954678E4, -1.32089234440210967447E6,
  80. 2.50083801823357915839E7, -5.38695755929454629881E7};
  81. static double PI180 = 1.74532925199432957692E-2;
  82. static double lossth = 1.0e14;
  83. #endif
  84. #ifdef DEC
  85. static unsigned short P[] = {0143514, 0113306, 0111171, 0174674,
  86. 0045214, 0147545, 0027744, 0167346,
  87. 0146210, 0177526, 0114514, 0105660};
  88. static unsigned short Q[] = {
  89. /*0040200,0000000,0000000,0000000,*/
  90. 0043525, 0142457, 0072633, 0025617, 0145241, 0036742, 0140525, 0162256,
  91. 0046276, 0146176, 0013526, 0143573, 0146515, 0077401, 0162762, 0150607};
  92. static unsigned short P1[] = {0036616, 0175065, 0011224, 0164711};
  93. #define PI180 *(double *)P1
  94. static double lossth = 8.0e14;
  95. #endif
  96. #ifdef IBMPC
  97. static unsigned short P[] = {0x3f38, 0xd24f, 0x92d8, 0xc0c9, 0x9ddd, 0xa5fc,
  98. 0x99ec, 0x4131, 0x9176, 0xd329, 0x1fea, 0xc171};
  99. static unsigned short Q[] = {
  100. /*0x0000,0x0000,0x0000,0x3ff0,*/
  101. 0x6572, 0xeeb3, 0xb8a5, 0x40ca, 0xbc96, 0x582a, 0x27bc, 0xc134,
  102. 0xd8ef, 0xc2ea, 0xd98f, 0x4177, 0x5a31, 0x3cbe, 0xafe0, 0xc189};
  103. static unsigned short P1[] = {0x9d39, 0xa252, 0xdf46, 0x3f91};
  104. #define PI180 *(double *)P1
  105. static double lossth = 1.0e14;
  106. #endif
  107. #ifdef MIEEE
  108. static unsigned short P[] = {0xc0c9, 0x92d8, 0xd24f, 0x3f38, 0x4131, 0x99ec,
  109. 0xa5fc, 0x9ddd, 0xc171, 0x1fea, 0xd329, 0x9176};
  110. static unsigned short Q[] = {0x40ca, 0xb8a5, 0xeeb3, 0x6572, 0xc134, 0x27bc,
  111. 0x582a, 0xbc96, 0x4177, 0xd98f, 0xc2ea, 0xd8ef,
  112. 0xc189, 0xafe0, 0x3cbe, 0x5a31};
  113. static unsigned short P1[] = {0x3f91, 0xdf46, 0xa252, 0x9d39};
  114. #define PI180 *(double *)P1
  115. static double lossth = 1.0e14;
  116. #endif
  117. #ifdef ANSIPROT
  118. extern double polevl(double, void *, int);
  119. extern double p1evl(double, void *, int);
  120. extern double floor(double);
  121. extern double ldexp(double, int);
  122. static double tancot(double, int);
  123. #else
  124. double polevl(), p1evl(), floor(), ldexp();
  125. static double tancot();
  126. #endif
  127. extern double MAXNUM;
  128. extern double PIO4;
  129. double tandg(x) double x;
  130. { return (tancot(x, 0)); }
  131. double cotdg(x) double x;
  132. { return (tancot(x, 1)); }
  133. static double tancot(xx, cotflg) double xx;
  134. int cotflg;
  135. {
  136. double x, y, z, zz;
  137. int j, sign;
  138. /* make argument positive but save the sign */
  139. if (xx < 0) {
  140. x = -xx;
  141. sign = -1;
  142. } else {
  143. x = xx;
  144. sign = 1;
  145. }
  146. if (x > lossth) {
  147. mtherr("tandg", TLOSS);
  148. return (0.0);
  149. }
  150. /* compute x mod PIO4 */
  151. y = floor(x / 45.0);
  152. /* strip high bits of integer part */
  153. z = ldexp(y, -3);
  154. z = floor(z); /* integer part of y/8 */
  155. z = y - ldexp(z, 3); /* y - 16 * (y/16) */
  156. /* integer and fractional part modulo one octant */
  157. j = z;
  158. /* map zeros and singularities to origin */
  159. if (j & 1) {
  160. j += 1;
  161. y += 1.0;
  162. }
  163. z = x - y * 45.0;
  164. z *= PI180;
  165. zz = z * z;
  166. if (zz > 1.0e-14)
  167. y = z + z * (zz * polevl(zz, P, 2) / p1evl(zz, Q, 4));
  168. else
  169. y = z;
  170. if (j & 2) {
  171. if (cotflg)
  172. y = -y;
  173. else {
  174. if (y != 0.0) {
  175. y = -1.0 / y;
  176. } else {
  177. mtherr("tandg", SING);
  178. y = MAXNUM;
  179. }
  180. }
  181. } else {
  182. if (cotflg) {
  183. if (y != 0.0)
  184. y = 1.0 / y;
  185. else {
  186. mtherr("cotdg", SING);
  187. y = MAXNUM;
  188. }
  189. }
  190. }
  191. if (sign < 0)
  192. y = -y;
  193. return (y);
  194. }