home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / ATAN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1987-06-20  |  4.6 KB  |  261 lines
  1. /*                            atan.c
  2.  *
  3.  *    Inverse circular tangent
  4.  *      (arctangent)
  5.  *
  6.  *
  7.  *
  8.  * SYNOPSIS:
  9.  *
  10.  * double x, y, atan();
  11.  *
  12.  * y = atan( 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   range      # trials      peak         rms
  31.  *    DEC       0, 10        9500       2.8e-17     8.6e-18
  32.  *    IEEE      0, 10        5000       2.4e-16     7.6e-17
  33.  *
  34.  */
  35. /*                            atan2()
  36.  *
  37.  *    Quadrant correct inverse circular tangent
  38.  *
  39.  *
  40.  *
  41.  * SYNOPSIS:
  42.  *
  43.  * double x, y, z, atan2();
  44.  *
  45.  * z = atan2( x, y );
  46.  *
  47.  *
  48.  *
  49.  * DESCRIPTION:
  50.  *
  51.  * Returns radian angle between 0 and +2pi whose tangent
  52.  * is y/x.
  53.  *
  54.  *
  55.  *
  56.  * ACCURACY:
  57.  *
  58.  * See atan.c.
  59.  *
  60.  */
  61.  
  62. /*                            atan.c */
  63.  
  64.  
  65. /*
  66. Cephes Math Library Release 2.0:  April, 1987
  67. Copyright 1984, 1987 by Stephen L. Moshier
  68. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  69. */
  70.  
  71.  
  72. #include "mconf.h"
  73.  
  74. #ifdef UNK
  75. static double P[] = {
  76. -8.40980878064499716001E-1,
  77. -8.83860837023772394279E0,
  78. -2.18476213081316705724E1,
  79. -1.48307050340438946993E1
  80. };
  81. static double Q[] = {
  82. /* 1.00000000000000000000E0,*/
  83.  1.54974124675307267552E1,
  84.  6.27906555762653017263E1,
  85.  9.22381329856214406485E1,
  86.  4.44921151021319438465E1,
  87. };
  88.  
  89. /* tan( 3*pi/8 ) */
  90. static double T3P8 = 2.41421356237309504880;
  91.  
  92. /* tan( pi/8 ) */
  93. static double TP8 = 0.41421356237309504880;
  94. #endif
  95.  
  96. #ifdef DEC
  97. static short P[] = {
  98. 0140127,0045205,0153731,0030027,
  99. 0141015,0065360,0116105,0025210,
  100. 0141256,0143755,0127056,0105716,
  101. 0141155,0045221,0056235,0072437
  102. };
  103. static short Q[] = {
  104. /*0040200,0000000,0000000,0000000,*/
  105. 0041167,0172546,0143212,0126224,
  106. 0041573,0024641,0116611,0153210,
  107. 0041670,0074754,0110422,0127624,
  108. 0041461,0173755,0002566,0014374
  109. };
  110.  
  111. /* tan( 3*pi/8 ) = 2.41421356237309504880 */
  112. static short T3P8A[] = {040432,0101171,0114774,0167462,};
  113. #define T3P8 *(double *)T3P8A
  114.  
  115. /* tan( pi/8 ) = 0.41421356237309504880 */
  116. static short TP8A[] = {037724,011714,0147747,074622,};
  117. #define TP8 *(double *)TP8A
  118. #endif
  119.  
  120. #ifdef IBMPC
  121. static short P[] = {
  122. 0x2603,0xbafb,0xe950,0xbfea,
  123. 0xa551,0x1388,0xad5e,0xc021,
  124. 0xd17a,0xb5c5,0xd8fd,0xc035,
  125. 0xaea4,0x2b93,0xa952,0xc02d
  126. };
  127. static short Q[] = {
  128. /*0x0000,0x0000,0x0000,0x3ff0,*/
  129. 0x5592,0xd8d1,0xfeac,0x402e,
  130. 0x3ad1,0x33b1,0x6534,0x404f,
  131. 0x55f2,0x9222,0x0f3d,0x4057,
  132. 0xc31f,0xa0ae,0x3efd,0x4046
  133. };
  134.  
  135. /* tan( 3*pi/8 ) = 2.41421356237309504880 */
  136. static short T3P8A[] = {0x9de6,0x333f,0x504f,0x4003};
  137. #define T3P8 *(double *)T3P8A
  138.  
  139. /* tan( pi/8 ) = 0.41421356237309504880 */
  140. static short TP8A[] = {0xef32,0x99fc,0x8279,0x3fda};
  141. #define TP8 *(double *)TP8A
  142. #endif
  143.  
  144. #ifdef MIEEE
  145. static short P[] = {
  146. 0xbfea,0xe950,0xbafb,0x2603,
  147. 0xc021,0xad5e,0x1388,0xa551,
  148. 0xc035,0xd8fd,0xb5c5,0xd17a,
  149. 0xc02d,0xa952,0x2b93,0xaea4
  150. };
  151. static short Q[] = {
  152. 0x402e,0xfeac,0xd8d1,0x5592,
  153. 0x404f,0x6534,0x33b1,0x3ad1,
  154. 0x4057,0x0f3d,0x9222,0x55f2,
  155. 0x4046,0x3efd,0xa0ae,0xc31f,
  156. };
  157.  
  158. /* tan( 3*pi/8 ) = 2.41421356237309504880 */
  159. static short T3P8A[] = {
  160. 0x4003,0x504f,0x333f,0x9de6
  161. };
  162. #define T3P8 *(double *)T3P8A
  163.  
  164. /* tan( pi/8 ) = 0.41421356237309504880 */
  165. static short TP8A[] = {
  166. 0x3fda,0x8279,0x99fc,0xef32
  167. };
  168. #define TP8 *(double *)TP8A
  169. #endif
  170.  
  171.  
  172. double atan(x)
  173. double x;
  174. {
  175. extern double PIO2, PIO4;
  176. double y,z,p,q;
  177. double polevl(), p1evl();
  178. short sign;
  179.  
  180. /* make argument positive and save the sign */
  181. sign = 1;
  182. if( x < 0.0 )
  183.     {
  184.     sign = -1;
  185.     x = -x;
  186.     }
  187.  
  188. /* range reduction */
  189. if( x > T3P8 )
  190.     {
  191.     y = PIO2;
  192.     x = -( 1.0/x );
  193.     }
  194.  
  195. else if( x > TP8 )
  196.     {
  197.     y = PIO4;
  198.     x = (x-1.0)/(x+1.0);
  199.     }
  200. else
  201.     y = 0.0;
  202.  
  203. /* rational form in x**2 */
  204. z = x * x;
  205. y = y + ( polevl( z, P, 3 ) / p1evl( z, Q, 4 ) ) * z * x + x;
  206.  
  207. if( sign < 0 )
  208.     y = -y;
  209.  
  210. return(y);
  211. }
  212.  
  213. /*                            atan2    */
  214.  
  215. extern double PI, PIO2;
  216.  
  217. double atan2( x, y )
  218. double x, y;
  219. {
  220. double z, w;
  221. short code;
  222. double atan();
  223.  
  224.  
  225. code = 0;
  226.  
  227. if( x < 0.0 )
  228.     code = 2;
  229. if( y < 0.0 )
  230.     code |= 1;
  231.  
  232. if( x == 0.0 )
  233.     {
  234.     if( code & 1 )
  235.         return( 3.0*PIO2 );
  236.     if( y == 0.0 )
  237.         return( 0.0 );
  238.     return( PIO2 );
  239.     }
  240.  
  241. if( y == 0.0 )
  242.     {
  243.     if( code & 2 )
  244.         return( PI );
  245.     return( 0.0 );
  246.     }
  247.  
  248.  
  249. switch( code )
  250.     {
  251.     case 0: w = 0.0; break;
  252.     case 1: w = 2.0 * PI; break;
  253.     case 2:
  254.     case 3: w = PI; break;
  255.     }
  256.  
  257. z = atan( y/x );
  258.  
  259. return( w + z );
  260. }
  261.