home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The World of Computer Software
/
World_Of_Computer_Software-02-385-Vol-1of3.iso
/
d
/
dec92.zip
/
1012088A
< prev
next >
Wrap
Text File
|
1992-10-13
|
4KB
|
117 lines
/* 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[_D0] & _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);
}
}