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