Listing 5 pow using the _Dtest macro to isolate the special cases

/* pow function */
#include "xmath.h"
double logt();
#define log(x) logt(x)

double (pow) (double x, double y) {     /* compute x^y */
   double yi, yx = y, z;
   int n, xexp, zexp, neg = 0, errx = _Dtest(x), inc;
#define shuge HUGE_EXP
#define dhuge (double)HUGE_EXP
#define ln2 0.69314718055994530942
   static const double rthalf = 0.70710678118654752440;
   const int erry = _Dtest(yx);
   long x0, y0 = ((unsigned long *) &rthalf)[_D0];
   _Dvar xx;

   xexp = 0;
   if (x < DBL_MIN) {
   /* shift up into normalized range */
      x *= 1 / DBL_EPSILON;
      xexp = -DBL_MANT_DIG;
   }
   xx._D = x;
   xexp += ((xx._W[_DO] &_DMASK) >> _DOFF) - _DBIAS;
   if (0 <= errx || 0 < erry) {
      /* x == 0, INF, NAN; y == INF, NAN */
      if (errx == NAN || erry == NAN)
         z = errx == NAN ? x : y, errx = NAN;
      else if (erry == INF)
         if (errx != INF)    /* 0^INF, finite^INF */
            errx = xexp <= 0 ? (y < 0 ? INF : 0) : (x == 1 || x == -1) ? 
               NAN : (y < 0 ? 0 : INF);
         else if (y == 0.0)
            return (1.0); /* x^0, x not a NaN */
         else if (errx == INF) {
            /* INF^finite */
            errx = y < 0.0 ? 0 : INF;
            neg = x < 0 && erry == 0 && ((int) y & 1);
         } else              /* 0^finite */
            errx = y < 0 0 ? INF : 0;
      if (errx == 0)
         return (0.0);
      else if (errx == INF) { /* return -INF or INF */
         errno = ERANGE;
         return x;
      } else {                /* return NaN */
         errno = EDOM;
         return (z);
      }
   }
   neg = 0;
   if (0.0 >= x) {
      x = -x, neg = (int) y & 1;
      if (y != (int) y) {    /* negative^fractional */
         y = (int) y;
         /* assume it's a small error */
         errno = EDOM;
      }
   }
   /* On some architectures, inc can be obtained
    * efficiently from "<" but this was the fastest
    * version on HP720 */
   xx._D = x;
   xexp -= inc = (unsigned long)
      ((x0 [_ = xx._W[_D0] & ~_DMASK |
         (_DBIAS << _DOFF)) - y0) >> 31;
   xx._W[_D0] = x0 + (inc << _DOFF);
   /* if xx < sqrt(1/2) shift up by factor 2 */
   x = xx._D;
   n = 0, yx = 0;
   if (y <= -dhuge)
      zexp = xexp < 0 ? shuge : xexp == 0 ? 0 : -shuge;
   else if (dhuge <= y)
      zexp = xexp < 0 ? -shuge : xexp == 0 ? 0 : shuge;
   else {
      /* y*log2(x) may be reasonable */
      double dexp;
      long zl;

      /* Form yx = y*xexp-zl carefully */
      yi = (2 << 16) / DBL_EPSILON;
      yi = y > 0 ? (y + yi) - yi : yi - (yi - y);
      dexp = xexp;
      yi = ((yi * dexp - (zl = y * dexp)) + (y - yi) *
         dexp) * ln2;
      zexp = zl <= -shuge ? -shuge : zl < shuge ?
         zl : shuge;
      if ((n = (int) y) > SAFE_EXP || n < -SAFE_EXP)
         n = 0;
   }
   {                           /* compute z = xfrac^n *
                           * 2^yx * 2^zexp */
      yx = yi + log(x) * (y - n);
      if (n < 0)
         n = -n;
      /* Z *= xfrac^n */
      z = (n & 1) ? x : 1.0;
      for (yi = x; n > 1; z *= yi)   /* scale by x^2^n */
         do
             yi *= yi;
         while (((n >>= 1) & 1) == 0);
      if (0 <= _Exp(&yx, zexp)) {     /* 2^yx*2^zexp */
         errno = ERANGE;   /* underflow or overflow */
         z = yx;
      } else
         /* either n==0 or there is no chance of
          * under/overflow */
         yi = yx;
         z = y < 0 ? yi / z : z * yi;
      /* don't invert and multiply */
      return (neg ? -z : z);
   }
}
/* End of File */