/* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z)
* on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1
*/
-static long double P[] = {
+static const long double P[] = {
8.3319510773868690346226E-4L,
4.9000050881978028599627E-1L,
1.7500123722550302671919E0L,
1.4000100839971580279335E0L,
};
-static long double Q[] = {
+static const long double Q[] = {
/* 1.0000000000000000000000E0L,*/
5.2500282295834889175431E0L,
8.4000598057587009834666E0L,
/* A[i] = 2^(-i/32), rounded to IEEE long double precision.
* If i is even, A[i] + B[i/2] gives additional accuracy.
*/
-static long double A[33] = {
+static const long double A[33] = {
1.0000000000000000000000E0L,
9.7857206208770013448287E-1L,
9.5760328069857364691013E-1L,
5.1094857432705833910408E-1L,
5.0000000000000000000000E-1L,
};
-static long double B[17] = {
+static const long double B[17] = {
0.0000000000000000000000E0L,
2.6176170809902549338711E-20L,
-1.0126791927256478897086E-20L,
/* 2^x = 1 + x P(x),
* on the interval -1/32 <= x <= 0
*/
-static long double R[] = {
+static const long double R[] = {
1.5089970579127659901157E-5L,
1.5402715328927013076125E-4L,
1.3333556028915671091390E-3L,
static const long double MAXLOGL = 1.1356523406294143949492E4L;
static const long double MINLOGL = -1.13994985314888605586758E4L;
static const long double LOGE2L = 6.9314718055994530941723E-1L;
-static volatile long double z;
-static long double w, W, Wa, Wb, ya, yb, u;
static const long double huge = 0x1p10000L;
/* XXX Prevent gcc from erroneously constant folding this. */
-static volatile long double twom10000 = 0x1p-10000L;
+static const volatile long double twom10000 = 0x1p-10000L;
static long double reducl(long double);
static long double powil(long double, int);
/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
int i, nflg, iyflg, yoddint;
long e;
+ volatile long double z=0;
+ long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
- if (y == 0.0L)
- return 1.0L;
+ if (y == 0.0)
+ return 1.0;
if (isnan(x))
return x;
if (isnan(y))
return y;
- if (y == 1.0L)
+ if (y == 1.0)
return x;
// FIXME: this is wrong, see pow special cases in c99 F.9.4.4
- if (!isfinite(y) && (x == -1.0L || x == 1.0L) )
+ if (!isfinite(y) && (x == -1.0 || x == 1.0) )
return y - y; /* +-1**inf is NaN */
- if (x == 1.0L)
- return 1.0L;
+ if (x == 1.0)
+ return 1.0;
if (y >= LDBL_MAX) {
- if (x > 1.0L)
+ if (x > 1.0)
return INFINITY;
- if (x > 0.0L && x < 1.0L)
- return 0.0L;
- if (x < -1.0L)
+ if (x > 0.0 && x < 1.0)
+ return 0.0;
+ if (x < -1.0)
return INFINITY;
- if (x > -1.0L && x < 0.0L)
- return 0.0L;
+ if (x > -1.0 && x < 0.0)
+ return 0.0;
}
if (y <= -LDBL_MAX) {
- if (x > 1.0L)
- return 0.0L;
- if (x > 0.0L && x < 1.0L)
+ if (x > 1.0)
+ return 0.0;
+ if (x > 0.0 && x < 1.0)
return INFINITY;
- if (x < -1.0L)
- return 0.0L;
- if (x > -1.0L && x < 0.0L)
+ if (x < -1.0)
+ return 0.0;
+ if (x > -1.0 && x < 0.0)
return INFINITY;
}
if (x >= LDBL_MAX) {
- if (y > 0.0L)
+ if (y > 0.0)
return INFINITY;
- return 0.0L;
+ return 0.0;
}
w = floorl(y);
yoddint = 0;
if (iyflg) {
ya = fabsl(y);
- ya = floorl(0.5L * ya);
- yb = 0.5L * fabsl(w);
+ ya = floorl(0.5 * ya);
+ yb = 0.5 * fabsl(w);
if( ya != yb )
yoddint = 1;
}
if (x <= -LDBL_MAX) {
- if (y > 0.0L) {
+ if (y > 0.0) {
if (yoddint)
return -INFINITY;
return INFINITY;
}
- if (y < 0.0L) {
+ if (y < 0.0) {
if (yoddint)
- return -0.0L;
+ return -0.0;
return 0.0;
}
}
nflg = 0; /* flag = 1 if x<0 raised to integer power */
- if (x <= 0.0L) {
- if (x == 0.0L) {
+ if (x <= 0.0) {
+ if (x == 0.0) {
if (y < 0.0) {
if (signbit(x) && yoddint)
return -INFINITY;
}
if (y > 0.0) {
if (signbit(x) && yoddint)
- return -0.0L;
+ return -0.0;
return 0.0;
}
- if (y == 0.0L)
- return 1.0L; /* 0**0 */
- return 0.0L; /* 0**y */
+ if (y == 0.0)
+ return 1.0; /* 0**0 */
+ return 0.0; /* 0**y */
}
if (iyflg == 0)
return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
*/
z = x*x;
w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3));
- w = w - ldexpl(z, -1); /* w - 0.5 * z */
+ w = w - 0.5*z;
/* Convert to base 2 logarithm:
* multiply by log2(e) = 1 + LOG2EA
/* Compute exponent term of the base 2 logarithm. */
w = -i;
- w = ldexpl(w, -LNXT); /* divide by NXT */
+ // TODO: use w * 0x1p-5;
+ w = scalbnl(w, -LNXT); /* divide by NXT */
w += e;
/* Now base 2 log of x is w + z. */
H = Fb + Gb;
Ha = reducl(H);
- w = ldexpl( Ga+Ha, LNXT );
+ w = scalbnl( Ga+Ha, LNXT );
/* Test the power of 2 for overflow */
if (w > MEXP)
e = w;
Hb = H - Ha;
- if (Hb > 0.0L) {
+ if (Hb > 0.0) {
e += 1;
- Hb -= 1.0L/NXT; /*0.0625L;*/
+ Hb -= 1.0/NXT; /*0.0625L;*/
}
/* Now the product y * log2(x) = Hb + e/NXT.
w = douba(e);
z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = z + w;
- z = ldexpl(z, i); /* multiply by integer power of 2 */
+ z = scalbnl(z, i); /* multiply by integer power of 2 */
if (nflg) {
/* For negative x,
* find out if the integer exponent
* is odd or even.
*/
- w = ldexpl(y, -1);
+ w = 0.5*y;
w = floorl(w);
- w = ldexpl(w, 1);
+ w = 2.0*w;
if (w != y)
z = -z; /* odd exponent */
}
{
long double t;
- t = ldexpl(x, LNXT);
+ t = scalbnl(x, LNXT);
t = floorl(t);
- t = ldexpl(t, -LNXT);
+ t = scalbnl(t, -LNXT);
return t;
}
long double s;
int n, e, sign, asign, lx;
- if (x == 0.0L) {
+ if (x == 0.0) {
if (nn == 0)
- return 1.0L;
+ return 1.0;
else if (nn < 0)
return LDBL_MAX;
- return 0.0L;
+ return 0.0;
}
if (nn == 0)
- return 1.0L;
+ return 1.0;
- if (x < 0.0L) {
+ if (x < 0.0) {
asign = -1;
x = -x;
} else
e = (lx - 1)*n;
if ((e == 0) || (e > 64) || (e < -64)) {
s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L);
- s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
+ s = (2.9142135623730950L * s - 0.5 + lx) * nn * LOGE2L;
} else {
s = LOGE2L * e;
}
* since roundoff error in 1.0/x will be amplified.
* The precise demarcation should be the gradual underflow threshold.
*/
- if (s < -MAXLOGL+2.0L) {
- x = 1.0L/x;
+ if (s < -MAXLOGL+2.0) {
+ x = 1.0/x;
sign = -sign;
}
if (n & 1)
y = x;
else {
- y = 1.0L;
+ y = 1.0;
asign = 0;
}
if (asign)
y = -y; /* odd power of negative number */
if (sign < 0)
- y = 1.0L/y;
+ y = 1.0/y;
return y;
}