lib/math/tan.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        <math.h>
  10 #include        <float.h>
  11 #include        <errno.h>
  12 #include        "localmath.h"
  13
  14 double
  15 tan(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         int negative = x < 0;
  23         int invert = 0;
  24         double  y;
  25         static double   p[] = {
  26                  1.0,
  27                 -0.13338350006421960681e+0,
  28                  0.34248878235890589960e-2,
  29                 -0.17861707342254426711e-4
  30         };
  31         static double   q[] = {
  32                  1.0,
  33                 -0.46671683339755294240e+0,
  34                  0.25663832289440112864e-1,
  35                 -0.31181531907010027307e-3,
  36                  0.49819433993786512270e-6
  37         };
  38
  39         if (__IsNan(x)) {
  40                 errno = EDOM;
  41                 return x;
  42         }
  43         if (negative) x = -x;
  44
  45         /* ??? avoid loss of significance, error if x is too large ??? */
  46
  47         y = x * M_2_PI + 0.5;
  48
  49         if (y >= DBL_MAX/M_PI_2) return 0.0;
  50
  51         /*      Use extended precision to calculate reduced argument.
  52                 Here we used 12 bits of the mantissa for a1.
  53                 Also split x in integer part x1 and fraction part x2.
  54         */
  55     #define A1 1.57080078125
  56     #define A2 -4.454455103380768678308e-6
  57         {
  58                 double x1, x2;
  59
  60                 modf(y, &y);
  61                 if (modf(0.5*y, &x1)) invert = 1;
  62                 x2 = modf(x, &x1);
  63                 x = x1 - y * A1;
  64                 x += x2;
  65                 x -= y * A2;
  66     #undef A1
  67     #undef A2
  68         }
  69
  70         /* ??? avoid underflow ??? */
  71         y = x * x;
  72         x += x * y * POLYNOM2(y, p+1);
  73         y = POLYNOM4(y, q);
  74         if (negative) x = -x;
  75         return invert ? -y/x : x/y;
  76 }