sys/src/libc/port/atan.c

   1 /*
   2         floating-point arctangent
   3
   4         atan returns the value of the arctangent of its
   5         argument in the range [-pi/2,pi/2].
   6
   7         atan2 returns the arctangent of arg1/arg2
   8         in the range [-pi,pi].
   9
  10         there are no error returns.
  11
  12         coefficients are #5077 from Hart & Cheney. (19.56D)
  13 */
  14
  15 #include <u.h>
  16 #include <libc.h>
  17
  18 #define sq2p1 2.414213562373095048802e0
  19 #define sq2m1  .414213562373095048802e0
  20 #define p4  .161536412982230228262e2
  21 #define p3  .26842548195503973794141e3
  22 #define p2  .11530293515404850115428136e4
  23 #define p1  .178040631643319697105464587e4
  24 #define p0  .89678597403663861959987488e3
  25 #define q4  .5895697050844462222791e2
  26 #define q3  .536265374031215315104235e3
  27 #define q2  .16667838148816337184521798e4
  28 #define q1  .207933497444540981287275926e4
  29 #define q0  .89678597403663861962481162e3
  30
  31
  32 /*
  33         xatan evaluates a series valid in the
  34         range [-0.414...,+0.414...]. (tan(pi/8))
  35  */
  36
  37 static
  38 double
  39 xatan(double arg)
  40 {
  41         double argsq, value;
  42
  43         argsq = arg*arg;
  44         value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0);
  45         value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0);
  46         return value*arg;
  47 }
  48
  49 /*
  50         satan reduces its argument (known to be positive)
  51         to the range [0,0.414...] and calls xatan.
  52  */
  53
  54 static
  55 double
  56 satan(double arg)
  57 {
  58
  59         if(arg < sq2m1)
  60                 return xatan(arg);
  61         if(arg > sq2p1)
  62                 return PIO2 - xatan(1/arg);
  63         return PIO2/2 + xatan((arg-1)/(arg+1));
  64 }
  65
  66 /*
  67         atan makes its argument positive and
  68         calls the inner routine satan.
  69  */
  70
  71 double
  72 atan(double arg)
  73 {
  74
  75         if(arg > 0)
  76                 return satan(arg);
  77         return -satan(-arg);
  78 }