Print this page
11210 libm should be cstyle(1ONBLD) clean
*** 20,29 ****
--- 20,30 ----
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
*** 38,52 ****
#define i1 1
#define i2 2
#define i3 3
static const long double zero = 0.0L, one = 1.0L, two = 2.0L;
-
extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
!
! static const long double
! two113 = 10384593717069655257060992658440192.0L,
ln2hi = 6.931471805599453094172319547495844850203e-0001L,
ln2lo = 1.667085920830552208890449330400379754169e-0025L,
A2 = 6.666666666666666666666666666666091393804e-0001L,
A3 = 4.000000000000000000000000407167070220671e-0001L,
A4 = 2.857142857142857142730077490612903681164e-0001L,
--- 39,50 ----
#define i1 1
#define i2 2
#define i3 3
static const long double zero = 0.0L, one = 1.0L, two = 2.0L;
extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
! static const long double two113 = 10384593717069655257060992658440192.0L,
ln2hi = 6.931471805599453094172319547495844850203e-0001L,
ln2lo = 1.667085920830552208890449330400379754169e-0025L,
A2 = 6.666666666666666666666666666666091393804e-0001L,
A3 = 4.000000000000000000000000407167070220671e-0001L,
A4 = 2.857142857142857142730077490612903681164e-0001L,
*** 62,126 ****
B7 = 1.333333333506731842033180638329317108428e-0001L,
B8 = 1.176469984587418890634302788283946761670e-0001L,
B9 = 1.053794891561452331722969901564862497132e-0001L;
static long double
! logl_x(long double x, long double *w) {
long double f, f1, v, s, z, qn, h, t;
! int *px = (int *) &x;
! int *pz = (int *) &z;
int i, j, ix, n;
n = 0;
ix = px[i0];
if (ix > 0x3ffef03f && ix < 0x3fff0820) { /* 65/63 > x > 63/65 */
f = x - one;
z = f * f;
if (((ix - 0x3fff0000) | px[i1] | px[i2] | px[i3]) == 0) {
*w = zero;
return (zero); /* log(1)= +0 */
}
qn = one / (two + f);
s = f * qn; /* |s|<2**-6 */
v = s * s;
! h = (long double) (2.0 * (double) s);
! f1 = (long double) ((double) f);
! t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn +
! s * (v * (B1 + v * (B2 + v * (B3 + v * (B4 +
! v * (B5 + v * (B6 + v * (B7 + v * (B8 + v * B9)))))))));
! s = (long double) ((double) (h + t));
*w = t - (s - h);
return (s);
}
if (ix < 0x00010000) { /* subnormal x */
x *= two113;
n = -113;
ix = px[i0];
}
/* LARGE_N */
n += ((ix + 0x200) >> 16) - 0x3fff;
ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */
px[i0] = ix;
i = ix + 0x200;
pz[i0] = i & 0xfffffc00;
pz[i1] = pz[i2] = pz[i3] = 0;
qn = one / (x + z);
f = x - z;
s = f * qn;
! f1 = (long double) ((double) f);
! h = (long double) (2.0 * (double) s);
t = qn * ((two * (f - z * h) - h * f1) - h * (f - f1));
j = (i >> 10) & 0x3f;
v = s * s;
! qn = (long double) n;
t += qn * ln2lo + _TBL_logl_lo[j];
! t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 +
! v * A7))))));
v = qn * ln2hi + _TBL_logl_hi[j];
s = h + v;
t += (h - (s - v));
! z = (long double) ((double) (s + t));
*w = t - (z - s);
return (z);
}
extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
--- 60,130 ----
B7 = 1.333333333506731842033180638329317108428e-0001L,
B8 = 1.176469984587418890634302788283946761670e-0001L,
B9 = 1.053794891561452331722969901564862497132e-0001L;
static long double
! logl_x(long double x, long double *w)
! {
long double f, f1, v, s, z, qn, h, t;
! int *px = (int *)&x;
! int *pz = (int *)&z;
int i, j, ix, n;
n = 0;
ix = px[i0];
+
if (ix > 0x3ffef03f && ix < 0x3fff0820) { /* 65/63 > x > 63/65 */
f = x - one;
z = f * f;
+
if (((ix - 0x3fff0000) | px[i1] | px[i2] | px[i3]) == 0) {
*w = zero;
return (zero); /* log(1)= +0 */
}
+
qn = one / (two + f);
s = f * qn; /* |s|<2**-6 */
v = s * s;
! h = (long double)(2.0 * (double)s);
! f1 = (long double)((double)f);
! t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn + s * (v *
! (B1 + v * (B2 + v * (B3 + v * (B4 + v * (B5 + v * (B6 + v *
! (B7 + v * (B8 + v * B9)))))))));
! s = (long double)((double)(h + t));
*w = t - (s - h);
return (s);
}
+
if (ix < 0x00010000) { /* subnormal x */
x *= two113;
n = -113;
ix = px[i0];
}
+
/* LARGE_N */
n += ((ix + 0x200) >> 16) - 0x3fff;
ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */
px[i0] = ix;
i = ix + 0x200;
pz[i0] = i & 0xfffffc00;
pz[i1] = pz[i2] = pz[i3] = 0;
qn = one / (x + z);
f = x - z;
s = f * qn;
! f1 = (long double)((double)f);
! h = (long double)(2.0 * (double)s);
t = qn * ((two * (f - z * h) - h * f1) - h * (f - f1));
j = (i >> 10) & 0x3f;
v = s * s;
! qn = (long double)n;
t += qn * ln2lo + _TBL_logl_lo[j];
! t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 + v *
! A7))))));
v = qn * ln2hi + _TBL_logl_hi[j];
s = h + v;
t += (h - (s - v));
! z = (long double)((double)(s + t));
*w = t - (z - s);
return (z);
}
extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
*** 129,192 ****
ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
ln2_32lo = 5.209643502595475652782654157501186731779e-27L,
ln2_64 = 1.083042469624914545964425189778400898568e-2L;
long double
! powl(long double x, long double y) {
long double z, ax;
long double y1, y2, w1, w2;
int sbx, sby, j, k, yisint, m;
int hx, lx, hy, ly, ahx, ahy;
! int *pz = (int *) &z;
! int *px = (int *) &x;
! int *py = (int *) &y;
hx = px[i0];
lx = px[i1] | px[i2] | px[i3];
hy = py[i0];
ly = py[i1] | py[i2] | py[i3];
ahx = hx & ~0x80000000;
ahy = hy & ~0x80000000;
if ((ahy | ly) == 0)
return (one); /* x**+-0 = 1 */
! else if (hx == 0x3fff0000 && lx == 0 &&
! (__xpg6 & _C99SUSv3_pow) != 0)
return (one); /* C99: 1**anything = 1 */
! else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) ||
! ahy > 0x7fff0000 || (ahy == 0x7fff0000 && ly != 0))
return (x + y); /* +-NaN return x+y */
/* includes Sun: 1**NaN = NaN */
! sbx = (unsigned) hx >> 31;
! sby = (unsigned) hy >> 31;
ax = fabsl(x);
/*
* 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 (sbx) {
! if (ahy >= 0x40700000) /* if |y|>=2**113 */
yisint = 2; /* even integer y */
! else if (ahy >= 0x3fff0000) {
k = (ahy >> 16) - 0x3fff; /* exponent */
if (k > 80) {
! j = ((unsigned) py[i3]) >> (112 - k);
if ((j << (112 - k)) == py[i3])
yisint = 2 - (j & 1);
} else if (k > 48) {
! j = ((unsigned) py[i2]) >> (80 - k);
if ((j << (80 - k)) == py[i2])
yisint = 2 - (j & 1);
} else if (k > 16) {
! j = ((unsigned) py[i1]) >> (48 - k);
if ((j << (48 - k)) == py[i1])
yisint = 2 - (j & 1);
} else if (ly == 0) {
j = ahy >> (16 - k);
if ((j << (16 - k)) == ahy)
yisint = 2 - (j & 1);
}
}
}
--- 133,204 ----
ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
ln2_32lo = 5.209643502595475652782654157501186731779e-27L,
ln2_64 = 1.083042469624914545964425189778400898568e-2L;
long double
! powl(long double x, long double y)
! {
long double z, ax;
long double y1, y2, w1, w2;
int sbx, sby, j, k, yisint, m;
int hx, lx, hy, ly, ahx, ahy;
! int *pz = (int *)&z;
! int *px = (int *)&x;
! int *py = (int *)&y;
hx = px[i0];
lx = px[i1] | px[i2] | px[i3];
hy = py[i0];
ly = py[i1] | py[i2] | py[i3];
ahx = hx & ~0x80000000;
ahy = hy & ~0x80000000;
if ((ahy | ly) == 0)
return (one); /* x**+-0 = 1 */
! else if (hx == 0x3fff0000 && lx == 0 && (__xpg6 & _C99SUSv3_pow) != 0)
return (one); /* C99: 1**anything = 1 */
! else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) || ahy >
! 0x7fff0000 || (ahy == 0x7fff0000 && ly != 0))
return (x + y); /* +-NaN return x+y */
+
/* includes Sun: 1**NaN = NaN */
! sbx = (unsigned)hx >> 31;
! sby = (unsigned)hy >> 31;
ax = fabsl(x);
+
/*
* 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 (sbx) {
! if (ahy >= 0x40700000) { /* if |y|>=2**113 */
yisint = 2; /* even integer y */
! } else if (ahy >= 0x3fff0000) {
k = (ahy >> 16) - 0x3fff; /* exponent */
+
if (k > 80) {
! j = ((unsigned)py[i3]) >> (112 - k);
!
if ((j << (112 - k)) == py[i3])
yisint = 2 - (j & 1);
} else if (k > 48) {
! j = ((unsigned)py[i2]) >> (80 - k);
!
if ((j << (80 - k)) == py[i2])
yisint = 2 - (j & 1);
} else if (k > 16) {
! j = ((unsigned)py[i1]) >> (48 - k);
!
if ((j << (48 - k)) == py[i1])
yisint = 2 - (j & 1);
} else if (ly == 0) {
j = ahy >> (16 - k);
+
if ((j << (16 - k)) == ahy)
yisint = 2 - (j & 1);
}
}
}
*** 198,221 ****
if ((__xpg6 & _C99SUSv3_pow) != 0)
return (one);
/* C99: (-1)**+-inf = 1 */
else
return (y - y);
/* Sun: (+-1)**+-inf = NaN */
! } else if (ahx >= 0x3fff0000)
/* (|x|>1)**+,-inf = inf,0 */
return (sby == 0 ? y : zero);
! else /* (|x|<1)**-,+inf = inf,0 */
return (sby != 0 ? -y : zero);
} else if (ahy == 0x3fff0000) { /* y is +-1 */
if (sby != 0)
return (one / x);
else
return (x);
! } else if (hy == 0x40000000) /* y is 2 */
return (x * x);
! else if (hy == 0x3ffe0000) { /* y is 0.5 */
if (!((ahx | lx) == 0 || ((ahx - 0x7fff0000) | lx) ==
0))
return (sqrtl(x));
}
}
--- 210,235 ----
if ((__xpg6 & _C99SUSv3_pow) != 0)
return (one);
/* C99: (-1)**+-inf = 1 */
else
return (y - y);
+
/* Sun: (+-1)**+-inf = NaN */
! } else if (ahx >= 0x3fff0000) {
/* (|x|>1)**+,-inf = inf,0 */
return (sby == 0 ? y : zero);
! } else { /* (|x|<1)**-,+inf = inf,0 */
return (sby != 0 ? -y : zero);
+ }
} else if (ahy == 0x3fff0000) { /* y is +-1 */
if (sby != 0)
return (one / x);
else
return (x);
! } else if (hy == 0x40000000) { /* y is 2 */
return (x * x);
! } else if (hy == 0x3ffe0000) { /* y is 0.5 */
if (!((ahx | lx) == 0 || ((ahx - 0x7fff0000) | lx) ==
0))
return (sqrtl(x));
}
}
*** 223,265 ****
/* special value of x */
if (lx == 0) {
if (ahx == 0x7fff0000 || ahx == 0 || ahx == 0x3fff0000) {
/* x is +-0,+-inf,+-1 */
z = ax;
if (sby == 1)
z = one / z; /* z = 1/|x| if y is negative */
if (sbx == 1) {
if (ahx == 0x3fff0000 && yisint == 0)
z = zero / zero;
/* (-1)**non-int is NaN */
else if (yisint == 1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
return (z);
}
}
/* (x<0)**(non-int) is NaN */
if (sbx == 1 && yisint == 0)
return (zero / zero); /* should be volatile */
! /* Now ax is finite, y is finite */
! /* first compute log(ax) = w1+w2, with 53 bits w1 */
w1 = logl_x(ax, &w2);
/* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
if (ly == 0 || ahy >= 0x43fe0000) {
y1 = y * w1;
y2 = y * w2;
} else {
! y1 = (long double) ((double) y);
y2 = (y - y1) * w1 + y * w2;
y1 *= w1;
}
z = y1 + y2;
j = pz[i0];
! if ((unsigned) j >= 0xffff0000) { /* NaN or -inf */
if (sbx == 1 && yisint == 1)
return (one / z);
else
return (-one / z);
} else if ((j & ~0x80000000) < 0x3fc30000) { /* |x|<2^-60 */
--- 237,286 ----
/* special value of x */
if (lx == 0) {
if (ahx == 0x7fff0000 || ahx == 0 || ahx == 0x3fff0000) {
/* x is +-0,+-inf,+-1 */
z = ax;
+
if (sby == 1)
z = one / z; /* z = 1/|x| if y is negative */
+
if (sbx == 1) {
if (ahx == 0x3fff0000 && yisint == 0)
z = zero / zero;
/* (-1)**non-int is NaN */
else if (yisint == 1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
+
return (z);
}
}
/* (x<0)**(non-int) is NaN */
if (sbx == 1 && yisint == 0)
return (zero / zero); /* should be volatile */
! /*
! * Now ax is finite, y is finite
! * first compute log(ax) = w1+w2, with 53 bits w1
! */
w1 = logl_x(ax, &w2);
/* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
if (ly == 0 || ahy >= 0x43fe0000) {
y1 = y * w1;
y2 = y * w2;
} else {
! y1 = (long double)((double)y);
y2 = (y - y1) * w1 + y * w2;
y1 *= w1;
}
+
z = y1 + y2;
j = pz[i0];
!
! if ((unsigned)j >= 0xffff0000) { /* NaN or -inf */
if (sbx == 1 && yisint == 1)
return (one / z);
else
return (-one / z);
} else if ((j & ~0x80000000) < 0x3fc30000) { /* |x|<2^-60 */
*** 272,318 ****
if (sbx == 1 && yisint == 1)
return (scalbnl(-one, 20000));
else
return (scalbnl(one, 20000));
}
! k = (int) (invln2_32 * (z + ln2_64));
} else {
! if ((unsigned) j > 0xc00d0000) {
if (sbx == 1 && yisint == 1)
return (scalbnl(-one, -20000));
else
return (scalbnl(one, -20000));
}
! k = (int) (invln2_32 * (z - ln2_64));
}
j = k & 0x1f;
m = k >> 5;
{
/* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
! long double
! t1 = 1.666666666666666666666666666660876387437e-1L,
t2 = -2.777777777777777777777707812093173478756e-3L,
t3 = 6.613756613756613482074280932874221202424e-5L,
t4 = -1.653439153392139954169609822742235851120e-6L,
t5 = 4.175314851769539751387852116610973796053e-8L;
! long double t = (long double) k;
w1 = (y2 - (t * ln2_32hi - y1)) - t * ln2_32lo;
t = w1 * w1;
w2 = (w1 - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) -
two;
z = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (w1 + w1)) / w2 -
_TBL_expl_lo[j]);
}
j = m + (pz[i0] >> 16);
! if (j && (unsigned) j < 0x7fff)
pz[i0] += m << 16;
else
z = scalbnl(z, m);
if (sbx == 1 && yisint == 1)
z = -z; /* (-ve)**(odd int) */
return (z);
}
#else
#error Unsupported Architecture
#endif /* defined(__sparc) */
--- 293,343 ----
if (sbx == 1 && yisint == 1)
return (scalbnl(-one, 20000));
else
return (scalbnl(one, 20000));
}
!
! k = (int)(invln2_32 * (z + ln2_64));
} else {
! if ((unsigned)j > 0xc00d0000) {
if (sbx == 1 && yisint == 1)
return (scalbnl(-one, -20000));
else
return (scalbnl(one, -20000));
}
!
! k = (int)(invln2_32 * (z - ln2_64));
}
+
j = k & 0x1f;
m = k >> 5;
{
/* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
! long double t1 = 1.666666666666666666666666666660876387437e-1L,
t2 = -2.777777777777777777777707812093173478756e-3L,
t3 = 6.613756613756613482074280932874221202424e-5L,
t4 = -1.653439153392139954169609822742235851120e-6L,
t5 = 4.175314851769539751387852116610973796053e-8L;
! long double t = (long double)k;
w1 = (y2 - (t * ln2_32hi - y1)) - t * ln2_32lo;
t = w1 * w1;
w2 = (w1 - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) -
two;
z = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (w1 + w1)) / w2 -
_TBL_expl_lo[j]);
}
j = m + (pz[i0] >> 16);
!
! if (j && (unsigned)j < 0x7fff)
pz[i0] += m << 16;
else
z = scalbnl(z, m);
if (sbx == 1 && yisint == 1)
z = -z; /* (-ve)**(odd int) */
+
return (z);
}
#else
#error Unsupported Architecture
#endif /* defined(__sparc) */