home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / TANL.C < prev    next >
C/C++ Source or Header  |  1993-01-30  |  5KB  |  252 lines
  1. /*                            tanl.c
  2.  *
  3.  *    Circular tangent, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, tanl();
  10.  *
  11.  * y = tanl( 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.  *    IEEE     +-1.07e9       30000     1.9e-19     4.8e-20
  30.  *
  31.  * ERROR MESSAGES:
  32.  *
  33.  *   message         condition          value returned
  34.  * tan total loss   x > 2^39                0.0
  35.  *
  36.  */
  37. /*                            cotl.c
  38.  *
  39.  *    Circular cotangent, long double precision
  40.  *
  41.  *
  42.  *
  43.  * SYNOPSIS:
  44.  *
  45.  * long double x, y, cotl();
  46.  *
  47.  * y = cotl( x );
  48.  *
  49.  *
  50.  *
  51.  * DESCRIPTION:
  52.  *
  53.  * Returns the circular cotangent of the radian argument x.
  54.  *
  55.  * Range reduction is modulo pi/4.  A rational function
  56.  *       x + x**3 P(x**2)/Q(x**2)
  57.  * is employed in the basic interval [0, pi/4].
  58.  *
  59.  *
  60.  *
  61.  * ACCURACY:
  62.  *
  63.  *                      Relative error:
  64.  * arithmetic   domain     # trials      peak         rms
  65.  *    IEEE     +-1.07e9      30000      1.9e-19     5.1e-20
  66.  *
  67.  *
  68.  * ERROR MESSAGES:
  69.  *
  70.  *   message         condition          value returned
  71.  * cot total loss   x > 2^39                0.0
  72.  * cot singularity  x = 0                  MAXNUM
  73.  *
  74.  */
  75.  
  76. /*
  77. Cephes Math Library Release 2.2:  December, 1990
  78. Copyright 1984, 1990 by Stephen L. Moshier
  79. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  80. */
  81.  
  82. #ifndef NOINCS
  83. #include "mconf.h"
  84. #endif
  85.  
  86. #ifdef UNK
  87. long double tanP[] = {
  88. -1.3093693918138377764608E4L,
  89.  1.1535166483858741613983E6L,
  90. -1.7956525197648487798769E7L,
  91. };
  92. long double tanQ[] = {
  93. /* 1.0000000000000000000000E0L,*/
  94.  1.3681296347069295467845E4L,
  95. -1.3208923444021096744731E6L,
  96.  2.5008380182335791583922E7L,
  97. -5.3869575592945462988123E7L,
  98. };
  99. long double tanDP1 = 7.853981554508209228515625E-1L;
  100. long double tanDP2 = 7.946627356147928367136046290398E-9L;
  101. long double tanDP3 = 3.061616997868382943065164830688E-17L;
  102. long double tanlossth = 5.49755813888e11L; /* 2^39 */
  103. #endif
  104.  
  105.  
  106. #ifdef IBMPC
  107. short tanP[] = {
  108. 0xbc1c,0x79f9,0xc692,0xcc96,0xc00c, XPD
  109. 0xe5b1,0xe4ee,0x652f,0x8ccf,0x4013, XPD
  110. 0xaf9a,0x4c8b,0x5699,0x88ff,0xc017, XPD
  111. };
  112. short tanQ[] = {
  113. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  114. 0x8ed4,0x9b2b,0x2f75,0xd5c5,0x400c, XPD
  115. 0xadcd,0x55e4,0xe2c1,0xa13d,0xc013, XPD
  116. 0x7adf,0x56c7,0x7e17,0xbecc,0x4017, XPD
  117. 0x86f6,0xf2d1,0x01e5,0xcd7f,0xc018, XPD
  118. };
  119. short tanP1[] = {0x0000,0x0000,0xda80,0xc90f,0x3ffe, XPD};
  120. short tanP2[] = {0x0000,0x0000,0xa300,0x8885,0x3fe4, XPD};
  121. short tanP3[] = {0x3707,0xa2e0,0x3198,0x8d31,0x3fc8, XPD};
  122. #define tanDP1 *(long double *)tanP1
  123. #define tanDP2 *(long double *)tanP2
  124. #define tanDP3 *(long double *)tanP3
  125. short tanlosi[] = {0x0000,0x0000,0x0000,0x8000,0x4026, XPD};
  126. #define tanlossth *(long double *)tanlosi
  127. #endif
  128.  
  129. #ifdef MIEEE
  130. long tanP[] = {
  131. 0xc00c0000,0xcc96c692,0x79f9bc1c,
  132. 0x40130000,0x8ccf652f,0xe4eee5b1,
  133. 0xc0170000,0x88ff5699,0x4c8baf9a,
  134. };
  135. long tanQ[] = {
  136. /*0x3fff0000,0x80000000,0x00000000,*/
  137. 0x400c0000,0xd5c52f75,0x9b2b8ed4,
  138. 0xc0130000,0xa13de2c1,0x55e4adcd,
  139. 0x40170000,0xbecc7e17,0x56c77adf,
  140. 0xc0180000,0xcd7f01e5,0xf2d186f6,
  141. };
  142. long tanP1[] = {0x3ffe0000,0xc90fda80,0x00000000};
  143. long tanP2[] = {0x3fe40000,0x8885a300,0x00000000};
  144. long tanP3[] = {0x3fc80000,0x8d313198,0xa2e03707};
  145. #define tanDP1 *(long double *)tanP1
  146. #define tanDP2 *(long double *)tanP2
  147. #define tanDP3 *(long double *)tanP3
  148. long tanlosi[] = {0x40260000,0x80000000,0x00000000};
  149. #define tanlossth *(long double *)tanlosi
  150. #endif
  151.  
  152. extern long double PIO4L, MACHEPL, MAXNUML;
  153. extern long double Zero, One;
  154.  
  155. long double tancotl();
  156.  
  157. long double tanl(x)
  158. long double x;
  159. {
  160.  
  161. return( tancotl(x,0) );
  162. }
  163.  
  164.  
  165. long double cotl(x)
  166. long double x;
  167. {
  168.  
  169. if( x == 0 )
  170.     {
  171.     mtherr( "cotl", SING );
  172.     return( MAXNUML );
  173.     }
  174. return( tancotl(x,1) );
  175. }
  176.  
  177.  
  178. long double tancotl( xx, cotflg )
  179. long double xx;
  180. int cotflg;
  181. {
  182. long double x, y, z, zz;
  183. int j, sign;
  184. long double polevll(), p1evll(), floorl(), ldexpl();
  185.  
  186. /* make argument positive but save the sign */
  187. if( xx < 0 )
  188.     {
  189.     x = -xx;
  190.     sign = -1;
  191.     }
  192. else
  193.     {
  194.     x = xx;
  195.     sign = 1;
  196.     }
  197.  
  198. if( x > tanlossth )
  199.     {
  200.     if( cotflg )
  201.         mtherr( "cotl", TLOSS );
  202.     else
  203.         mtherr( "tanl", TLOSS );
  204.     return(Zero);
  205.     }
  206.  
  207. /* compute x mod PIO4 */
  208. y = floorl( x/PIO4L );
  209.  
  210. /* strip high bits of integer part */
  211. z = ldexpl( y, -4 );
  212. z = floorl(z);        /* integer part of y/16 */
  213. z = y - ldexpl( z, 4 );  /* y - 16 * (y/16) */
  214.  
  215. /* integer and fractional part modulo one octant */
  216. j = (int )z;
  217.  
  218. /* map zeros and singularities to origin */
  219. if( j & 1 )
  220.     {
  221.     j += 1;
  222.     y += One;
  223.     }
  224.  
  225. z = ((x - y * tanDP1) - y * tanDP2) - y * tanDP3;
  226.  
  227. zz = z * z;
  228.  
  229. if( zz > MACHEPL )
  230.     y = z  +  z * (zz * polevll( zz, tanP, 2 )/p1evll(zz, tanQ, 4));
  231. else
  232.     y = z;
  233.     
  234. if( j & 2 )
  235.     {
  236.     if( cotflg )
  237.         y = -y;
  238.     else
  239.         y = -One/y;
  240.     }
  241. else
  242.     {
  243.     if( cotflg )
  244.         y = One/y;
  245.     }
  246.  
  247. if( sign < 0 )
  248.     y = -y;
  249.  
  250. return( y );
  251. }
  252.