home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cmath / tan.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  5.1 KB  |  283 lines
  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      0,+1.07e9      44000     4.0e-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                  MAXNUM
  74.  *
  75.  */
  76.  
  77. /*
  78. Cephes Math Library Release 2.1:  March, 1989
  79. Copyright 1984, 1987, 1989 by Stephen L. Moshier
  80. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  81. */
  82.  
  83. #include "mconf.h"
  84.  
  85. #ifdef UNK
  86. static double P[] = {
  87. -1.30936939181383777646E4,
  88.  1.15351664838587416140E6,
  89. -1.79565251976484877988E7
  90. };
  91. static double Q[] = {
  92. /* 1.00000000000000000000E0,*/
  93.  1.36812963470692954678E4,
  94. -1.32089234440210967447E6,
  95.  2.50083801823357915839E7,
  96. -5.38695755929454629881E7
  97. };
  98. static double DP1 = 7.853981554508209228515625E-1;
  99. static double DP2 = 7.94662735614792836714E-9;
  100. static double DP3 = 3.06161699786838294307E-17;
  101. static double lossth = 1.073741824e9;
  102. #endif
  103.  
  104. #ifdef DEC
  105. static short P[] = {
  106. 0143514,0113306,0111171,0174674,
  107. 0045214,0147545,0027744,0167346,
  108. 0146210,0177526,0114514,0105660
  109. };
  110. static short Q[] = {
  111. /*0040200,0000000,0000000,0000000,*/
  112. 0043525,0142457,0072633,0025617,
  113. 0145241,0036742,0140525,0162256,
  114. 0046276,0146176,0013526,0143573,
  115. 0146515,0077401,0162762,0150607
  116. };
  117. /*  7.853981629014015197753906250000E-1 */
  118. static short P1[] = {0040111,0007732,0120000,0000000,};
  119. /*  4.960467869796758577649598009884E-10 */
  120. static short P2[] = {0030410,0055060,0100000,0000000,};
  121. /*  2.860594363054915898381331279295E-18 */
  122. static short P3[] = {0021523,0011431,0105056,0001560,};
  123. #define DP1 *(double *)P1
  124. #define DP2 *(double *)P2
  125. #define DP3 *(double *)P3
  126. static double lossth = 1.073741824e9;
  127. #endif
  128.  
  129. #ifdef IBMPC
  130. static short P[] = {
  131. 0x3f38,0xd24f,0x92d8,0xc0c9,
  132. 0x9ddd,0xa5fc,0x99ec,0x4131,
  133. 0x9176,0xd329,0x1fea,0xc171
  134. };
  135. static short Q[] = {
  136. /*0x0000,0x0000,0x0000,0x3ff0,*/
  137. 0x6572,0xeeb3,0xb8a5,0x40ca,
  138. 0xbc96,0x582a,0x27bc,0xc134,
  139. 0xd8ef,0xc2ea,0xd98f,0x4177,
  140. 0x5a31,0x3cbe,0xafe0,0xc189
  141. };
  142. /*
  143.   7.85398125648498535156E-1,
  144.   3.77489470793079817668E-8,
  145.   2.69515142907905952645E-15,
  146. */
  147. static short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
  148. static short P2[] = {0x0000,0x0000,0x442d,0x3e64};
  149. static short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
  150. #define DP1 *(double *)P1
  151. #define DP2 *(double *)P2
  152. #define DP3 *(double *)P3
  153. static double lossth = 1.073741824e9;
  154. #endif
  155.  
  156. #ifdef MIEEE
  157. static short P[] = {
  158. 0xc0c9,0x92d8,0xd24f,0x3f38,
  159. 0x4131,0x99ec,0xa5fc,0x9ddd,
  160. 0xc171,0x1fea,0xd329,0x9176
  161. };
  162. static short Q[] = {
  163. 0x40ca,0xb8a5,0xeeb3,0x6572,
  164. 0xc134,0x27bc,0x582a,0xbc96,
  165. 0x4177,0xd98f,0xc2ea,0xd8ef,
  166. 0xc189,0xafe0,0x3cbe,0x5a31
  167. };
  168. static short P1[] = {
  169. 0x3fe9,0x21fb,0x4000,0x0000
  170. };
  171. static short P2[] = {
  172. 0x3e64,0x442d,0x0000,0x0000
  173. };
  174. static short P3[] = {
  175. 0x3ce8,0x4698,0x98cc,0x5170,
  176. };
  177. #define DP1 *(double *)P1
  178. #define DP2 *(double *)P2
  179. #define DP3 *(double *)P3
  180. static double lossth = 1.073741824e9;
  181. #endif
  182.  
  183. extern double MAXNUM;
  184. extern double PIO4;
  185.  
  186. double tan(x)
  187. double x;
  188. {
  189. double tancot();
  190.  
  191. return( tancot(x,0) );
  192. }
  193.  
  194.  
  195. double cot(x)
  196. double x;
  197. {
  198. double tancot();
  199.  
  200. if( x == 0.0 )
  201.     {
  202.     mtherr( "cot", SING );
  203.     return( MAXNUM );
  204.     }
  205. return( tancot(x,1) );
  206. }
  207.  
  208.  
  209. static double tancot( xx, cotflg )
  210. double xx;
  211. int cotflg;
  212. {
  213. double x, y, z, zz;
  214. int j, sign;
  215. double polevl(), p1evl(), floor(), ldexp();
  216.  
  217. /* make argument positive but save the sign */
  218. if( xx < 0 )
  219.     {
  220.     x = -xx;
  221.     sign = -1;
  222.     }
  223. else
  224.     {
  225.     x = xx;
  226.     sign = 1;
  227.     }
  228.  
  229. if( x > lossth )
  230.     {
  231.     if( cotflg )
  232.         mtherr( "cot", TLOSS );
  233.     else
  234.         mtherr( "tan", TLOSS );
  235.     return(0.0);
  236.     }
  237.  
  238. /* compute x mod PIO4 */
  239. y = floor( x/PIO4 );
  240.  
  241. /* strip high bits of integer part */
  242. z = ldexp( y, -3 );
  243. z = floor(z);        /* integer part of y/8 */
  244. z = y - ldexp( z, 3 );  /* y - 16 * (y/16) */
  245.  
  246. /* integer and fractional part modulo one octant */
  247. j = z;
  248.  
  249. /* map zeros and singularities to origin */
  250. if( j & 1 )
  251.     {
  252.     j += 1;
  253.     y += 1.0;
  254.     }
  255.  
  256. z = ((x - y * DP1) - y * DP2) - y * DP3;
  257.  
  258. zz = z * z;
  259.  
  260. if( zz > 1.0e-14 )
  261.     y = z  +  z * (zz * polevl( zz, P, 2 )/p1evl(zz, Q, 4));
  262. else
  263.     y = z;
  264.     
  265. if( j & 2 )
  266.     {
  267.     if( cotflg )
  268.         y = -y;
  269.     else
  270.         y = -1.0/y;
  271.     }
  272. else
  273.     {
  274.     if( cotflg )
  275.         y = 1.0/y;
  276.     }
  277.  
  278. if( sign < 0 )
  279.     y = -y;
  280.  
  281. return( y );
  282. }
  283.