home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ldouble / tanl.c < prev   
Encoding:
C/C++ Source or Header  |  1992-11-17  |  4.6 KB  |  247 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. #include "mconf.h"
  83.  
  84. #ifdef UNK
  85. static long double P[] = {
  86. -1.3093693918138377764608E4L,
  87.  1.1535166483858741613983E6L,
  88. -1.7956525197648487798769E7L,
  89. };
  90. static long double Q[] = {
  91. /* 1.0000000000000000000000E0L,*/
  92.  1.3681296347069295467845E4L,
  93. -1.3208923444021096744731E6L,
  94.  2.5008380182335791583922E7L,
  95. -5.3869575592945462988123E7L,
  96. };
  97. static long double DP1 = 7.853981554508209228515625E-1L;
  98. static long double DP2 = 7.946627356147928367136046290398E-9L;
  99. static long double DP3 = 3.061616997868382943065164830688E-17L;
  100. #endif
  101.  
  102.  
  103. #ifdef IBMPC
  104. static short P[] = {
  105. 0xbc1c,0x79f9,0xc692,0xcc96,0xc00c,
  106. 0xe5b1,0xe4ee,0x652f,0x8ccf,0x4013,
  107. 0xaf9a,0x4c8b,0x5699,0x88ff,0xc017,
  108. };
  109. static short Q[] = {
  110. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  111. 0x8ed4,0x9b2b,0x2f75,0xd5c5,0x400c,
  112. 0xadcd,0x55e4,0xe2c1,0xa13d,0xc013,
  113. 0x7adf,0x56c7,0x7e17,0xbecc,0x4017,
  114. 0x86f6,0xf2d1,0x01e5,0xcd7f,0xc018,
  115. };
  116. static short P1[] = {0x0000,0x0000,0xda80,0xc90f,0x3ffe};
  117. static short P2[] = {0x0000,0x0000,0xa300,0x8885,0x3fe4};
  118. static short P3[] = {0x3707,0xa2e0,0x3198,0x8d31,0x3fc8};
  119. #define DP1 *(long double *)P1
  120. #define DP2 *(long double *)P2
  121. #define DP3 *(long double *)P3
  122. #endif
  123.  
  124. #ifdef MIEEE
  125. static long P[] = {
  126. 0xc00c0000,0xcc96c692,0x79f9bc1c,
  127. 0x40130000,0x8ccf652f,0xe4eee5b1,
  128. 0xc0170000,0x88ff5699,0x4c8baf9a,
  129. };
  130. static long Q[] = {
  131. /*0x3fff0000,0x80000000,0x00000000,*/
  132. 0x400c0000,0xd5c52f75,0x9b2b8ed4,
  133. 0xc0130000,0xa13de2c1,0x55e4adcd,
  134. 0x40170000,0xbecc7e17,0x56c77adf,
  135. 0xc0180000,0xcd7f01e5,0xf2d186f6,
  136. };
  137. static long P1[] = {0x3ffe0000,0xc90fda80,0x00000000};
  138. static long P2[] = {0x3fe40000,0x8885a300,0x00000000};
  139. static long P3[] = {0x3fc80000,0x8d313198,0xa2e03707};
  140. #define DP1 *(long double *)P1
  141. #define DP2 *(long double *)P2
  142. #define DP3 *(long double *)P3
  143. #endif
  144.  
  145. static long double lossth = 5.49755813888e11L; /* 2^39 */
  146. extern long double PIO4L;
  147. extern long double MAXNUML;
  148.  
  149.  
  150. long double tanl(x)
  151. long double x;
  152. {
  153. long double tancotl();
  154.  
  155. return( tancotl(x,0) );
  156. }
  157.  
  158.  
  159. long double cotl(x)
  160. long double x;
  161. {
  162. long double tancotl();
  163.  
  164. if( x == 0.0L )
  165.     {
  166.     mtherr( "cotl", SING );
  167.     return( MAXNUML );
  168.     }
  169. return( tancotl(x,1) );
  170. }
  171.  
  172.  
  173. static long double tancotl( xx, cotflg )
  174. long double xx;
  175. int cotflg;
  176. {
  177. long double x, y, z, zz;
  178. int j, sign;
  179. long double polevll(), p1evll(), floorl(), ldexpl();
  180.  
  181. /* make argument positive but save the sign */
  182. if( xx < 0.0L )
  183.     {
  184.     x = -xx;
  185.     sign = -1;
  186.     }
  187. else
  188.     {
  189.     x = xx;
  190.     sign = 1;
  191.     }
  192.  
  193. if( x > lossth )
  194.     {
  195.     if( cotflg )
  196.         mtherr( "cotl", TLOSS );
  197.     else
  198.         mtherr( "tanl", TLOSS );
  199.     return(0.0L);
  200.     }
  201.  
  202. /* compute x mod PIO4 */
  203. y = floorl( x/PIO4L );
  204.  
  205. /* strip high bits of integer part */
  206. z = ldexpl( y, -4 );
  207. z = floorl(z);        /* integer part of y/16 */
  208. z = y - ldexpl( z, 4 );  /* y - 16 * (y/16) */
  209.  
  210. /* integer and fractional part modulo one octant */
  211. j = z;
  212.  
  213. /* map zeros and singularities to origin */
  214. if( j & 1 )
  215.     {
  216.     j += 1;
  217.     y += 1.0L;
  218.     }
  219.  
  220. z = ((x - y * DP1) - y * DP2) - y * DP3;
  221.  
  222. zz = z * z;
  223.  
  224. if( zz > 1.0e-20L )
  225.     y = z  +  z * (zz * polevll( zz, P, 2 )/p1evll(zz, Q, 4));
  226. else
  227.     y = z;
  228.     
  229. if( j & 2 )
  230.     {
  231.     if( cotflg )
  232.         y = -y;
  233.     else
  234.         y = -1.0L/y;
  235.     }
  236. else
  237.     {
  238.     if( cotflg )
  239.         y = 1.0L/y;
  240.     }
  241.  
  242. if( sign < 0 )
  243.     y = -y;
  244.  
  245. return( y );
  246. }
  247.