lib/math/atan.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7 /* $Header$ */
   8
   9 #include        <float.h>
  10 #include        <math.h>
  11 #include        <errno.h>
  12 #include        "localmath.h"
  13
  14 double
  15 atan(double x)
  16 {
  17         /*      Algorithm and coefficients from:
  18                         "Software manual for the elementary functions"
  19                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  20         */
  21
  22         static double p[] = {
  23                 -0.13688768894191926929e+2,
  24                 -0.20505855195861651981e+2,
  25                 -0.84946240351320683534e+1,
  26                 -0.83758299368150059274e+0
  27         };
  28         static double q[] = {
  29                  0.41066306682575781263e+2,
  30                  0.86157349597130242515e+2,
  31                  0.59578436142597344465e+2,
  32                  0.15024001160028576121e+2,
  33                  1.0
  34         };
  35         static double a[] = {
  36                 0.0,
  37                 0.52359877559829887307710723554658381,  /* pi/6 */
  38                 M_PI_2,
  39                 1.04719755119659774615421446109316763   /* pi/3 */
  40         };
  41
  42         int     neg = x < 0;
  43         int     n;
  44         double  g;
  45
  46         if (__IsNan(x)) {
  47                 errno = EDOM;
  48                 return x;
  49         }
  50         if (neg) {
  51                 x = -x;
  52         }
  53         if (x > 1.0) {
  54                 x = 1.0/x;
  55                 n = 2;
  56         }
  57         else    n = 0;
  58
  59         if (x > 0.26794919243112270647) {       /* 2-sqtr(3) */
  60                 n = n + 1;
  61                 x = (((0.73205080756887729353*x-0.5)-0.5)+x)/
  62                         (1.73205080756887729353+x);
  63         }
  64
  65         /* ??? avoid underflow ??? */
  66
  67         g = x * x;
  68         x += x * g * POLYNOM3(g, p) / POLYNOM4(g, q);
  69         if (n > 1) x = -x;
  70         x += a[n];
  71         return neg ? -x : x;
  72 }