home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / ATANL.C < prev    next >
C/C++ Source or Header  |  1993-01-30  |  5KB  |  265 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. #ifndef NOINCS
  77. #include "mconf.h"
  78. #endif
  79.  
  80. #ifdef UNK
  81. long double atanP[] = {
  82. -8.6863818178092187535440E-1L,
  83. -1.4683508633175792446076E1L,
  84. -6.3976888655834347413154E1L,
  85. -9.9988763777265819915721E1L,
  86. -5.0894116899623603312185E1L,
  87. };
  88. long double atanQ[] = {
  89. /* 1.00000000000000000000E0L,*/
  90.  2.2981886733594175366172E1L,
  91.  1.4399096122250781605352E2L,
  92.  3.6144079386152023162701E2L,
  93.  3.9157570175111990631099E2L,
  94.  1.5268235069887081006606E2L,
  95. };
  96.  
  97. /* tan( 3*pi/8 ) */
  98. long double T3P8 =  2.41421356237309504880169L;
  99.  
  100. /* tan( pi/8 ) */
  101. long double TP8 =  4.1421356237309504880169L;
  102. #endif
  103.  
  104.  
  105. #ifdef IBMPC
  106. short atanP[] = {
  107. 0x8ece,0xce53,0x1266,0xde5f,0xbffe, XPD
  108. 0x07e6,0xa061,0xa6bf,0xeaef,0xc002, XPD
  109. 0x53ee,0xf291,0x557f,0xffe8,0xc004, XPD
  110. 0xf9d6,0xeda6,0x3f3e,0xc7fa,0xc005, XPD
  111. 0xb6c3,0x6abc,0x9361,0xcb93,0xc004, XPD
  112. };
  113. short atanQ[] = {
  114. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  115. 0x54d4,0x894e,0xe76e,0xb7da,0x4003, XPD
  116. 0x76b9,0x7a46,0xafa2,0x8ffd,0x4006, XPD
  117. 0xe3a9,0xe9c0,0x6bee,0xb4b8,0x4007, XPD
  118. 0xabc1,0x50a7,0xb098,0xc3c9,0x4007, XPD
  119. 0x891c,0x100d,0xae89,0x98ae,0x4006, XPD
  120. };
  121.  
  122. /* tan( 3*pi/8 ) = 2.41421356237309504880 */
  123. short T3P8A[] = {0x3242,0xfcef,0x7999,0x9a82,0x4000, XPD};
  124. #define T3P8 *(long double *)T3P8A
  125.  
  126. /* tan( pi/8 ) = 0.41421356237309504880 */
  127. short TP8A[] = {0x9211,0xe779,0xcccf,0xd413,0x3ffd, XPD};
  128. #define TP8 *(long double *)TP8A
  129. #endif
  130.  
  131. #ifdef MIEEE
  132. long atanP[] = {
  133. 0xbffe0000,0xde5f1266,0xce538ece,
  134. 0xc0020000,0xeaefa6bf,0xa06107e6,
  135. 0xc0040000,0xffe8557f,0xf29153ee,
  136. 0xc0050000,0xc7fa3f3e,0xeda6f9d6,
  137. 0xc0040000,0xcb939361,0x6abcb6c3,
  138. };
  139. long atanQ[] = {
  140. /*0x3fff0000,0x80000000,0x00000000,*/
  141. 0x40030000,0xb7dae76e,0x894e54d4,
  142. 0x40060000,0x8ffdafa2,0x7a4676b9,
  143. 0x40070000,0xb4b86bee,0xe9c0e3a9,
  144. 0x40070000,0xc3c9b098,0x50a7abc1,
  145. 0x40060000,0x98aeae89,0x100d891c,
  146. };
  147.  
  148. /* tan( 3*pi/8 ) = 2.41421356237309504880 */
  149. long T3P8A[] = {0x40000000,0x9a827999,0xfcef3242};
  150. #define T3P8 *(long double *)T3P8A
  151.  
  152. /* tan( pi/8 ) = 0.41421356237309504880 */
  153. long TP8A[] = {0x3ffd0000,0xd413cccf,0xe7799211};
  154. #define TP8 *(long double *)TP8A
  155. #endif
  156. extern long double Zero, One, Two, Three;
  157.  
  158. long double atanl(x)
  159. long double x;
  160. {
  161. extern long double PIO2L, PIO4L;
  162. long double y, z;
  163. long double polevll(), p1evll();
  164. short sign;
  165.  
  166. /* make argument positive and save the sign */
  167. sign = 1;
  168. if( x < 0 )
  169.     {
  170.     sign = -1;
  171.     x = -x;
  172.     }
  173.  
  174. /* range reduction */
  175. if( x > T3P8 )
  176.     {
  177.     y = PIO2L;
  178.     x = -( One/x );
  179.     }
  180.  
  181. else if( x > TP8 )
  182.     {
  183.     y = PIO4L;
  184.     x = (x-One)/(x+One);
  185.     }
  186. else
  187.     y = Zero;
  188.  
  189. /* rational form in x**2 */
  190. z = x * x;
  191. y = y + ( polevll( z, atanP, 4 ) / p1evll( z, atanQ, 5 ) ) * z * x + x;
  192.  
  193. if( sign < 0 )
  194.     y = -y;
  195.  
  196. return(y);
  197. }
  198.  
  199. /*                            atan2    */
  200.  
  201.  
  202. extern long double PIL, PIO2L;
  203.  
  204. #if ANSIC
  205. long double atan2l( y, x )
  206. #else
  207. long double atan2l( x, y )
  208. #endif
  209. long double x, y;
  210. {
  211. long double z, w;
  212. short code;
  213. long double atanl();
  214.  
  215.  
  216. code = 0;
  217.  
  218. if( x < 0 )
  219.     code = 2;
  220. if( y < 0 )
  221.     code |= 1;
  222.  
  223. if( x == 0 )
  224.     {
  225.     if( code & 1 )
  226.         {
  227. #if ANSIC
  228.         return( -PIO2L );
  229. #else
  230.         return( Three*PIO2L );
  231. #endif
  232.         }
  233.     if( y == 0 )
  234.         return( Zero );
  235.     return( PIO2L );
  236.     }
  237.  
  238. if( y == 0 )
  239.     {
  240.     if( code & 2 )
  241.         return( PIL );
  242.     return( Zero );
  243.     }
  244.  
  245. #if ANSIC
  246. if( code < 2 )
  247.     w = Zero;
  248. else if( code == 2 )
  249.     w = PIL;
  250. else
  251.     w = -PIL;
  252. #else
  253. if( code == 0 )
  254.     w = Zero;
  255. else if( code == 1 )
  256.     w = Two * PIL;
  257. else
  258.     w = PIL;
  259. #endif
  260.  
  261. z = atanl( y/x );
  262.  
  263. return( w + z );
  264. }
  265.