home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ldouble / atanl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  4.7 KB  |  264 lines
  1. /*                            atanl.c
  2.  *
  3.  *    Inverse circular tangent, long double precision
  4.  *      (arctangent)
  5.  *
  6.  *
  7.  *
  8.  * SYNOPSIS:
  9.  *
  10.  * long double x, y, atanl();
  11.  *
  12.  * y = atanl( x );
  13.  *
  14.  *
  15.  *
  16.  * DESCRIPTION:
  17.  *
  18.  * Returns radian angle between -pi/2 and +pi/2 whose tangent
  19.  * is x.
  20.  *
  21.  * Range reduction is from four intervals into the interval
  22.  * from zero to  tan( pi/8 ).  The approximant uses a rational
  23.  * function of degree 3/4 of the form x + x**3 P(x)/Q(x).
  24.  *
  25.  *
  26.  *
  27.  * ACCURACY:
  28.  *
  29.  *                      Relative error:
  30.  * arithmetic   domain     # trials      peak         rms
  31.  *    IEEE      -10, 10    150000       1.3e-19     3.0e-20
  32.  *
  33.  */
  34. /*                            atan2l()
  35.  *
  36.  *    Quadrant correct inverse circular tangent,
  37.  *    long double precision
  38.  *
  39.  *
  40.  *
  41.  * SYNOPSIS:
  42.  *
  43.  * long double x, y, z, atan2l();
  44.  *
  45.  * z = atan2l( y, x );
  46.  *
  47.  *
  48.  *
  49.  * DESCRIPTION:
  50.  *
  51.  * Returns radian angle whose tangent is y/x.
  52.  * Define compile time symbol ANSIC = 1 for ANSI standard,
  53.  * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
  54.  * 0 to 2PI, args (x,y).
  55.  *
  56.  *
  57.  *
  58.  * ACCURACY:
  59.  *
  60.  *                      Relative error:
  61.  * arithmetic   domain     # trials      peak         rms
  62.  *    IEEE      -10, 10     60000       1.7e-19     3.2e-20
  63.  * See atan.c.
  64.  *
  65.  */
  66.  
  67. /*                            atan.c */
  68.  
  69.  
  70. /*
  71. Cephes Math Library Release 2.2:  December, 1990
  72. Copyright 1984, 1990 by Stephen L. Moshier
  73. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  74. */
  75.  
  76.  
  77. #include "mconf.h"
  78.  
  79. #ifdef UNK
  80. static long double P[] = {
  81. -8.6863818178092187535440E-1L,
  82. -1.4683508633175792446076E1L,
  83. -6.3976888655834347413154E1L,
  84. -9.9988763777265819915721E1L,
  85. -5.0894116899623603312185E1L,
  86. };
  87. static long double Q[] = {
  88. /* 1.00000000000000000000E0L,*/
  89.  2.2981886733594175366172E1L,
  90.  1.4399096122250781605352E2L,
  91.  3.6144079386152023162701E2L,
  92.  3.9157570175111990631099E2L,
  93.  1.5268235069887081006606E2L,
  94. };
  95.  
  96. /* tan( 3*pi/8 ) */
  97. static long double T3P8 =  2.41421356237309504880169L;
  98.  
  99. /* tan( pi/8 ) */
  100. static long double TP8 =  4.1421356237309504880169e-1L;
  101. #endif
  102.  
  103.  
  104. #ifdef IBMPC
  105. static short P[] = {
  106. 0x8ece,0xce53,0x1266,0xde5f,0xbffe,
  107. 0x07e6,0xa061,0xa6bf,0xeaef,0xc002,
  108. 0x53ee,0xf291,0x557f,0xffe8,0xc004,
  109. 0xf9d6,0xeda6,0x3f3e,0xc7fa,0xc005,
  110. 0xb6c3,0x6abc,0x9361,0xcb93,0xc004,
  111. };
  112. static short Q[] = {
  113. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  114. 0x54d4,0x894e,0xe76e,0xb7da,0x4003,
  115. 0x76b9,0x7a46,0xafa2,0x8ffd,0x4006,
  116. 0xe3a9,0xe9c0,0x6bee,0xb4b8,0x4007,
  117. 0xabc1,0x50a7,0xb098,0xc3c9,0x4007,
  118. 0x891c,0x100d,0xae89,0x98ae,0x4006,
  119. };
  120.  
  121. /* tan( 3*pi/8 ) = 2.41421356237309504880 */
  122. static short T3P8A[] = {0x3242,0xfcef,0x7999,0x9a82,0x4000};
  123. #define T3P8 *(long double *)T3P8A
  124.  
  125. /* tan( pi/8 ) = 0.41421356237309504880 */
  126. static short TP8A[] = {0x9211,0xe779,0xcccf,0xd413,0x3ffd};
  127. #define TP8 *(long double *)TP8A
  128. #endif
  129.  
  130. #ifdef MIEEE
  131. static long P[] = {
  132. 0xbffe0000,0xde5f1266,0xce538ece,
  133. 0xc0020000,0xeaefa6bf,0xa06107e6,
  134. 0xc0040000,0xffe8557f,0xf29153ee,
  135. 0xc0050000,0xc7fa3f3e,0xeda6f9d6,
  136. 0xc0040000,0xcb939361,0x6abcb6c3,
  137. };
  138. static long Q[] = {
  139. /*0x3fff0000,0x80000000,0x00000000,*/
  140. 0x40030000,0xb7dae76e,0x894e54d4,
  141. 0x40060000,0x8ffdafa2,0x7a4676b9,
  142. 0x40070000,0xb4b86bee,0xe9c0e3a9,
  143. 0x40070000,0xc3c9b098,0x50a7abc1,
  144. 0x40060000,0x98aeae89,0x100d891c,
  145. };
  146.  
  147. /* tan( 3*pi/8 ) = 2.41421356237309504880 */
  148. static long T3P8A[] = {0x40000000,0x9a827999,fcef3242};
  149. #define T3P8 *(long double *)T3P8A
  150.  
  151. /* tan( pi/8 ) = 0.41421356237309504880 */
  152. static long TP8A[] = {0x3ffd0000,0xd413cccf,0xe7799211};
  153. #define TP8 *(long double *)TP8A
  154. #endif
  155.  
  156.  
  157. long double atanl(x)
  158. long double x;
  159. {
  160. extern long double PIO2L, PIO4L;
  161. long double y,z,p,q;
  162. long double polevll(), p1evll();
  163. short sign;
  164.  
  165. /* make argument positive and save the sign */
  166. sign = 1;
  167. if( x < 0.0L )
  168.     {
  169.     sign = -1;
  170.     x = -x;
  171.     }
  172.  
  173. /* range reduction */
  174. if( x > T3P8 )
  175.     {
  176.     y = PIO2L;
  177.     x = -( 1.0L/x );
  178.     }
  179.  
  180. else if( x > TP8 )
  181.     {
  182.     y = PIO4L;
  183.     x = (x-1.0L)/(x+1.0L);
  184.     }
  185. else
  186.     y = 0.0L;
  187.  
  188. /* rational form in x**2 */
  189. z = x * x;
  190. y = y + ( polevll( z, P, 4 ) / p1evll( z, Q, 5 ) ) * z * x + x;
  191.  
  192. if( sign < 0 )
  193.     y = -y;
  194.  
  195. return(y);
  196. }
  197.  
  198. /*                            atan2    */
  199.  
  200.  
  201. extern long double PIL, PIO2L;
  202.  
  203. #if ANSIC
  204. long double atan2l( y, x )
  205. #else
  206. long double atan2l( x, y )
  207. #endif
  208. long double x, y;
  209. {
  210. long double z, w;
  211. short code;
  212. long double atanl();
  213.  
  214.  
  215. code = 0;
  216.  
  217. if( x < 0.0L )
  218.     code = 2;
  219. if( y < 0.0L )
  220.     code |= 1;
  221.  
  222. if( x == 0.0L )
  223.     {
  224.     if( code & 1 )
  225.         {
  226. #if ANSIC
  227.         return( -PIO2L );
  228. #else
  229.         return( 3.0L*PIO2L );
  230. #endif
  231.         }
  232.     if( y == 0.0L )
  233.         return( 0.0L );
  234.     return( PIO2L );
  235.     }
  236.  
  237. if( y == 0.0L )
  238.     {
  239.     if( code & 2 )
  240.         return( PIL );
  241.     return( 0.0L );
  242.     }
  243.  
  244.  
  245. switch( code )
  246.     {
  247. #if ANSIC
  248.     case 0:
  249.     case 1: w = 0.0L; break;
  250.     case 2: w = PIL; break;
  251.     case 3: w = -PIL; break;
  252. #else
  253.     case 0: w = 0.0L; break;
  254.     case 1: w = 2.0L * PIL; break;
  255.     case 2:
  256.     case 3: w = PIL; break;
  257. #endif
  258.     }
  259.  
  260. z = atanl( y/x );
  261.  
  262. return( w + z );
  263. }
  264.