lib/ack/libp/atn.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
   8 /* $Header$ */
   9
  10 #define __NO_DEFS
  11 #include <math.h>
  12
  13 #if __STDC__
  14 #include <pc_math.h>
  15 #endif
  16
  17 double
  18 _atn(x)
  19         double x;
  20 {
  21         /*      Algorithm and coefficients from:
  22                         "Software manual for the elementary functions"
  23                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  24         */
  25
  26         static double p[] = {
  27                 -0.13688768894191926929e+2,
  28                 -0.20505855195861651981e+2,
  29                 -0.84946240351320683534e+1,
  30                 -0.83758299368150059274e+0
  31         };
  32         static double q[] = {
  33                  0.41066306682575781263e+2,
  34                  0.86157349597130242515e+2,
  35                  0.59578436142597344465e+2,
  36                  0.15024001160028576121e+2,
  37                  1.0
  38         };
  39         static double a[] = {
  40                 0.0,
  41                 0.52359877559829887307710723554658381,  /* pi/6 */
  42                 M_PI_2,
  43                 1.04719755119659774615421446109316763   /* pi/3 */
  44         };
  45
  46         int     neg = x < 0;
  47         int     n;
  48         double  g;
  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 }