/* @(#)e_pow.c 1.5 04/04/22 SMI */ /* * ==================================================== * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. * * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ====================================================
*/
//#include <sys/cdefs.h> //__FBSDID("$FreeBSD$");
/* __ieee754_pow(x,y) return x**y * * n * Method: Let x = 2 * (1+f) * 1. Compute and return log2(x) in two pieces: * log2(x) = w1 + w2, * where w1 has 53-24 = 29 bit trailing zeros. * 2. Perform y*log2(x) = n+y' by simulating multi-precision * arithmetic, where |y'|<=0.5. * 3. Return x**y = 2**n*exp(y'*log2) * * Special cases: * 1. (anything) ** 0 is 1 * 2. (anything) ** 1 is itself * 3. (anything) ** NAN is NAN except 1 ** NAN = 1 * 4. NAN ** (anything except 0) is NAN * 5. +-(|x| > 1) ** +INF is +INF * 6. +-(|x| > 1) ** -INF is +0 * 7. +-(|x| < 1) ** +INF is +0 * 8. +-(|x| < 1) ** -INF is +INF * 9. +-1 ** +-INF is 1 * 10. +0 ** (+anything except 0, NAN) is +0 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 * 12. +0 ** (-anything except 0, NAN) is +INF * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) * 15. +INF ** (+anything except 0,NAN) is +INF * 16. +INF ** (-anything except 0,NAN) is +0 * 17. -INF ** (anything) = -0 ** (-anything) * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) * 19. (-anything except 0 and inf) ** (non-integer) is NAN * * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular * pow(integer,integer) * always returns the correct integer provided it is * representable. * * Constants : * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown.
*/
/* x==1: 1**y = 1, even if y is NaN */ if (hx==0x3ff00000 && lx == 0) return one;
/* y!=zero: result is NaN if either arg is NaN */ if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) return nan_mix(x, y);
/* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer * yisint = 1 ... y is an odd int * yisint = 2 ... y is an even int
*/
yisint = 0; if(hx<0) { if(iy>=0x43400000) yisint = 2; /* even integer y */ elseif(iy>=0x3ff00000) {
k = (iy>>20)-0x3ff; /* exponent */ if(k>20) {
j = ly>>(52-k); if(((u_int32_t)j<<(52-k))==ly) yisint = 2-(j&1);
} elseif(ly==0) {
j = iy>>(20-k); if((j<<(20-k))==iy) yisint = 2-(j&1);
}
}
}
/* special value of y */ if(ly==0) { if (iy==0x7ff00000) { /* y is +-inf */ if(((ix-0x3ff00000)|lx)==0) return one; /* (-1)**+-inf is 1 */ elseif (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ return (hy>=0)? y: zero; else/* (|x|<1)**-,+inf = inf,0 */ return (hy<0)?-y: zero;
} if(iy==0x3ff00000) { /* y is +-1 */ if(hy<0) return one/x; elsereturn x;
} if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3fe00000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ return std::sqrt(x);
}
}
ax = fabs(x); /* special value of x */ if(lx==0) { if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
z = ax; /*x is +-0,+-inf,+-1*/ if(hy<0) z = one/z; /* z = (1/|x|) */ if(hx<0) { if(((ix-0x3ff00000)|yisint)==0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
} elseif(yisint==1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
} return z;
}
}
/* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be n = (hx>>31)+1; but ANSI C says a right shift of a signed negative quantity is
implementation defined. */
n = ((u_int32_t)hx>>31)-1;
/* (x<0)**(non-int) is NaN */ if((n|yisint)==0) return (x-x)/(x-x);
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
/* |y| is huge */ if(iy>0x41e00000) { /* if |y| > 2**31 */ if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
} /* over/underflow if x is not close to one */ if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny; if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny; /* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */
t = ax-one; /* t has 20 trailing zeros */
w = (t*t)*(half-t*(thrd-t*qrtr));
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;
t1 = u+v;
SET_LOW_WORD(t1,0);
t2 = v-(t1-u);
} else { double ss,s2,s_h,s_l,t_h,t_l;
n = 0; /* take care subnormal number */ if(ix<0x00100000)
{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff; /* determine interval */
ix = j|0x3ff00000; /* normalize ix */ if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ elseif(j<0xBB67A) k=1; /* |x|<sqrt(3) */ else {k=0;n+=1;ix -= 0x00100000;}
SET_HIGH_WORD(ax,ix);
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
