sys/src/ape/lib/ap/math/tan.c

   1 /*
   2         floating point tangent
   3
   4         A series is used after range reduction.
   5         Coefficients are #4285 from Hart & Cheney. (19.74D)
   6  */
   7
   8 #include <math.h>
   9 #include <errno.h>
  10
  11 static double invpi       = 1.27323954473516268;
  12 static double p0         = -0.1306820264754825668269611177e+5;
  13 static double p1          = 0.1055970901714953193602353981e+4;
  14 static double p2         = -0.1550685653483266376941705728e+2;
  15 static double p3          = 0.3422554387241003435328470489e-1;
  16 static double p4          = 0.3386638642677172096076369e-4;
  17 static double q0         = -0.1663895238947119001851464661e+5;
  18 static double q1          = 0.4765751362916483698926655581e+4;
  19 static double q2         = -0.1555033164031709966900124574e+3;
  20
  21 double
  22 tan(double arg)
  23 {
  24         double sign, temp, e, x, xsq;
  25         int flag, i;
  26
  27         flag = 0;
  28         sign = 1;
  29         if(arg < 0){
  30                 arg = -arg;
  31                 sign = -1;
  32         }
  33         arg = arg*invpi;   /* overflow? */
  34         x = modf(arg, &e);
  35         i = e;
  36         switch(i%4) {
  37         case 1:
  38                 x = 1 - x;
  39                 flag = 1;
  40                 break;
  41
  42         case 2:
  43                 sign = - sign;
  44                 flag = 1;
  45                 break;
  46
  47         case 3:
  48                 x = 1 - x;
  49                 sign = - sign;
  50                 break;
  51
  52         case 0:
  53                 break;
  54         }
  55
  56         xsq = x*x;
  57         temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x;
  58         temp = temp/(((xsq+q2)*xsq+q1)*xsq+q0);
  59
  60         if(flag == 1) {
  61                 if(temp == 0) {
  62                         errno = EDOM;
  63                         if (sign > 0)
  64                                 return HUGE_VAL;
  65                         return -HUGE_VAL;
  66                 }
  67                 temp = 1/temp;
  68         }
  69         return sign*temp;
  70 }