home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_10_12 / 1012088a < prev    next >
Text File  |  1992-10-13  |  4KB  |  117 lines
  1. /* pow function */
  2. #include "xmath.h"
  3. double logt();
  4. #define log(x) logt(x)
  5.  
  6. double (pow) (double x, double y) {     /* compute x^y */
  7.     double yi, yx = y, z;
  8.     int n, xexp, zexp, neg = 0, errx = _Dtest(x), inc;
  9. #define shuge HUGE_EXP
  10. #define dhuge (double)HUGE_EXP
  11. #define ln2 0.69314718055994530942
  12.     static const double rthalf = 0.70710678118654752440;
  13.     const int erry = _Dtest(yx);
  14.     long x0, y0 = ((unsigned long *) &rthalf)[_D0];
  15.     _Dvar xx;
  16.  
  17.     xexp = 0;
  18.     if (x < DBL_MIN) {
  19.     /* shift up into normalized range */
  20.         x *= 1 / DBL_EPSILON;
  21.         xexp = -DBL_MANT_DIG;
  22.     }
  23.     xx._D = x;
  24.     xexp += ((xx._W[_D0] & _DMASK) >> _DOFF) - _DBIAS;
  25.     if (0 <= errx || 0 < erry) {
  26.         /* x == 0, INF, NAN; y == INF, NAN */
  27.         if (errx == NAN || erry == NAN)
  28.             z = errx == NAN ? x : y, errx = NAN;
  29.         else if (erry == INF)
  30.             if (errx != INF)    /* 0^INF, finite^INF */
  31.                 errx = xexp <= 0 ? (y < 0 ? INF : 0)
  32.                     :
  33.                     (x == 1 || x == -1) ? NAN
  34.                     : (y < 0 ? 0 : INF);
  35.             else if (y == 0.0)
  36.                 return (1.0);   /* x^0, x not a NaN */
  37.             else if (errx == INF) {
  38.                 /* INF^finite */
  39.                 errx = y < 0.0 ? 0 : INF;
  40.                 neg = x < 0 && erry == 0 && ((int) y & 1);
  41.             } else              /* 0^finite */
  42.                 errx = y < 0.0 ? INF : 0;
  43.         if (errx == 0)
  44.             return (0.0);
  45.         else if (errx == INF) { /* return -INF or INF */
  46.             errno = ERANGE;
  47.             return x;
  48.         } else {                /* return NaN */
  49.             errno = EDOM;
  50.             return (z);
  51.         }
  52.     }
  53.     neg = 0;
  54.     if (0.0 >= x) {
  55.         x = -x, neg = (int) y & 1;
  56.         if (y != (int) y) {     /* negative^fractional */
  57.             y = (int) y;
  58.             /* assume it's a small error */
  59.             errno = EDOM;
  60.         }
  61.     }
  62.     /* On some architectures, inc can be obtained
  63.      * efficiently from "<" but this was the fastest
  64.      * version on HP720 */
  65.     xx._D = x;
  66.     xexp -= inc = (unsigned long)
  67.         ((x0 = xx._W[_D0] & ~_DMASK |
  68.             (_DBIAS << _DOFF)) - y0) >> 31;
  69.     xx._W[_D0] = x0 + (inc << _DOFF);
  70.     /* if xx < sqrt(1/2) shift up by factor 2 */
  71.     x = xx._D;
  72.     n = 0, yx = 0;
  73.     if (y <= -dhuge)
  74.         zexp = xexp < 0 ? shuge : xexp == 0 ? 0 : -shuge;
  75.     else if (dhuge <= y)
  76.         zexp = xexp < 0 ? -shuge : xexp == 0 ? 0 : shuge;
  77.     else {
  78.         /* y*log2(x) may be reasonable */
  79.         double dexp;
  80.         long zl;
  81.  
  82.         /* Form yx = y*xexp-zl carefully */
  83.         yi = (2 << 16) / DBL_EPSILON;
  84.         yi = y > 0 ? (y + yi) - yi : yi - (yi - y);
  85.         dexp = xexp;
  86.         yi = ((yi * dexp - (zl = y * dexp)) + (y - yi) *
  87.             dexp) * ln2;
  88.         zexp = zl <= -shuge ? -shuge : zl < shuge ?
  89.             zl : shuge;
  90.         if ((n = (int) y) > SAFE_EXP || n < -SAFE_EXP)
  91.             n = 0;
  92.     }
  93.     {                           /* compute z = xfrac^n *
  94.                                  * 2^yx * 2^zexp */
  95.         yx = yi + log(x) * (y - n);
  96.         if (n < 0)
  97.             n = -n;
  98.         /* Z *= xfrac^n */
  99.         z = (n & 1) ? x : 1.0;
  100.         for (yi = x; n > 1; z *= yi)    /* scale by x^2^n */
  101.             do
  102.                 yi *= yi;
  103.             while (((n >>= 1) & 1) == 0);
  104.         if (0 <= _Exp(&yx, zexp)) {     /* 2^yx*2^zexp */
  105.             errno = ERANGE;     /* underflow or overflow */
  106.             z = yx;
  107.         } else
  108.             /* either n==0 or there is no chance of
  109.              * under/overflow */
  110.             yi = yx;
  111.             z = y < 0 ? yi / z : z * yi;
  112.         /* don't invert and multiply */
  113.         return (neg ? -z : z);
  114.     }
  115. }
  116.  
  117.