Listing 6 (tan_2.c)

typedef struct {
   double X1, X2;
}      ARG_D_2;                 /* vector 2 */
ARG_D_2
tan_2(xi1, xi2)
   double xi1, xi2;
{
   double x2, x, n1, x4;
   ARG_D_2 res;
#include "float.h"
#if FLT_ROUNDS != 1
       #error "rounding mode not nearest; adjust code"
#endif
#if FLT_RADIX !=2 && FLT_RADIX != 10
       #error "code not optimum for accuracy in this
RADIX"
#endif
#include <errno.h>
#ifndef errno
   extern int errno;
#endif
#include <math.h>
#define M_PI  3.14159265358979323846
   x2 = (n1 = (xi1 > 0 ? 1 / LDBL_EPSILON :
             -1 / LDBL_EPSILON))7 + (x = xi1) / M_PI;
   x -= (x2 - n1) * M_PI;
   if (fabs(xi1 / M_PI) >= 1 / LDBL_EPSILON |
      fabs(xi2 / M_PI) >= 1 / LDBL_EPSILON)
      errno = ERANGE;
 /* now in 1st or 4th quadrant */
#define c0 33281881.3202530279
   n1 = c0 + (x2 = x * x) * (-15666569.8711211851);
   x4 = x2 * x2;
   res.X1 = x * (c0 + x2 * (-4572609.43103684572) + x4 *
     (131095.887915363619 + x2 * (-968.863245687503149 + x2))) / (n1 + x4 * (915701.668921990803
                       + x2 * (-13491.7937027796916)
                       + x4 * 44.4083322286368691));
 /* copy 2 */
   x2 = (n1 = (xi2 > 0 ? 1 / LDBL_EPSILON :
              -1 / LDBL_EPSILON)) + (x = xi2) / M_PI;
   x -= (x2 - n1) * M_PI;
   n1 = 915701.668921990803 - (x2 = x * x)
      * 13491.7937027796916;
   x4 = x2 * x2;
   res.X2 = (x * (c0 + x2 * (-4572609.43103684572)) +
     (131095.887915363619 + x2 * (-968.863245687503149 +
    x2)) * x4 * x) / (c0 + x2 * (-15666569.8711211851) +
                x4 * (n1 + x4 * 44.4083322286368691));
   return res;
}
/* End of File */