-/* @(#)e_pow.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
-#endif
-
/* __ieee754_pow(x,y) return x**y
*
* n
* 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 muti-precision
+ * 2. Perform y*log2(x) = n+y' by simulating muti-precision
* arithmetic, where |y'|<=0.5.
* 3. Return x**y = 2**n*exp(y'*log2)
*
* Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular
* pow(integer,integer)
- * always returns the correct integer provided it is
+ * 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
+ * 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.
*/
#include "math.h"
#include "math_private.h"
-#ifdef __STDC__
-static const double
-#else
-static double
-#endif
+static const double
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
-#ifdef __STDC__
- double __ieee754_pow(double x, double y)
-#else
- double __ieee754_pow(x,y)
- double x, y;
-#endif
+double attribute_hidden __ieee754_pow(double x, double y)
{
double z,ax,z_h,z_l,p_h,p_l;
double y1,t1,t2,r,s,t,u,v,w;
ix = hx&0x7fffffff; iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
- if((iy|ly)==0) return one;
+ if((iy|ly)==0) return one;
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
- iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
- return x+y;
+ iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
+ return x+y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 2 ... y is an even int
*/
yisint = 0;
- if(hx<0) {
+ if(hx<0) {
if(iy>=0x43400000) yisint = 2; /* even integer y */
else if(iy>=0x3ff00000) {
k = (iy>>20)-0x3ff; /* exponent */
j = iy>>(20-k);
if((j<<(20-k))==iy) yisint = 2-(j&1);
}
- }
- }
+ }
+ }
/* special value of y */
- if(ly==0) {
+ if(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if(((ix-0x3ff00000)|lx)==0)
return y - y; /* inf**+-1 is NaN */
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; else return x;
}
if(hy==0x40000000) return x*x; /* y is 2 */
if(hy==0x3fe00000) { /* y is 0.5 */
if(hx>=0) /* x >= +0 */
- return __ieee754_sqrt(x);
+ return __ieee754_sqrt(x);
}
}
if(hx<0) {
if(((ix-0x3ff00000)|yisint)==0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
- } else if(yisint==1)
+ } else if(yisint==1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
return z;
}
}
-
+
/* (x<0)**(non-int) is NaN */
if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
/* over/underflow if x is not close to one */
if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
- /* now |1-x| is tiny <= 2**-20, suffice to compute
+ /* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */
t = x-1; /* t has 20 trailing zeros */
w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
- }
+ }
t = p_l+p_h;
SET_LOW_WORD(t,0);
u = t*lg2_h;
else SET_HIGH_WORD(z,j);
return s*z;
}
+
+/*
+ * wrapper pow(x,y) return x**y
+ */
+#ifndef _IEEE_LIBM
+double pow(double x, double y)
+{
+ double z = __ieee754_pow(x, y);
+ if (_LIB_VERSION == _IEEE_|| isnan(y))
+ return z;
+ if (isnan(x)) {
+ if (y == 0.0)
+ return __kernel_standard(x, y, 42); /* pow(NaN,0.0) */
+ return z;
+ }
+ if (x == 0.0) {
+ if (y == 0.0)
+ return __kernel_standard(x, y, 20); /* pow(0.0,0.0) */
+ if (isfinite(y) && y < 0.0)
+ return __kernel_standard(x,y,23); /* pow(0.0,negative) */
+ return z;
+ }
+ if (!isfinite(z)) {
+ if (isfinite(x) && isfinite(y)) {
+ if (isnan(z))
+ return __kernel_standard(x, y, 24); /* pow neg**non-int */
+ return __kernel_standard(x, y, 21); /* pow overflow */
+ }
+ }
+ if (z == 0.0 && isfinite(x) && isfinite(y))
+ return __kernel_standard(x, y, 22); /* pow underflow */
+ return z;
+}
+#else
+strong_alias(__ieee754_pow, pow)
+#endif
+libm_hidden_def(pow)
OSDN Git Service
