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.
*/
*** 42,52 ****
* unsigned short exp:15;
* unsigned short frac1:16;
*/
#ifdef __LITTLE_ENDIAN
-
/* array indices used to access words within a double */
#define HIWORD 1
#define LOWORD 0
/* structure used to access words within a quad */
--- 43,52 ----
*** 59,79 ****
} l;
long double d;
};
/* default NaN returned for sqrt(neg) */
! static const union longdouble
! qnan = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff };
/* signalling NaN used to raise invalid */
static const union {
unsigned u[2];
double d;
} snan = { 0, 0x7ff00001 };
-
#else
-
/* array indices used to access words within a double */
#define HIWORD 0
#define LOWORD 1
/* structure used to access words within a quad */
--- 59,78 ----
} l;
long double d;
};
/* default NaN returned for sqrt(neg) */
! static const union longdouble qnan = {
! 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff
! };
/* signalling NaN used to raise invalid */
static const union {
unsigned u[2];
double d;
} snan = { 0, 0x7ff00001 };
#else
/* array indices used to access words within a double */
#define HIWORD 0
#define LOWORD 1
/* structure used to access words within a quad */
*** 86,109 ****
} l;
long double d;
};
/* default NaN returned for sqrt(neg) */
! static const union longdouble
! qnan = { 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff };
/* signalling NaN used to raise invalid */
static const union {
unsigned u[2];
double d;
} snan = { 0x7ff00001, 0 };
-
#endif /* __LITTLE_ENDIAN */
!
! static const double
! zero = 0.0,
half = 0.5,
one = 1.0,
huge = 1.0e300,
tiny = 1.0e-300,
two36 = 6.87194767360000000000e+10,
--- 85,106 ----
} l;
long double d;
};
/* default NaN returned for sqrt(neg) */
! static const union longdouble qnan = {
! 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff
! };
/* signalling NaN used to raise invalid */
static const union {
unsigned u[2];
double d;
} snan = { 0x7ff00001, 0 };
#endif /* __LITTLE_ENDIAN */
! static const double zero = 0.0,
half = 0.5,
one = 1.0,
huge = 1.0e300,
tiny = 1.0e-300,
two36 = 6.87194767360000000000e+10,
*** 118,185 ****
twom66 = 1.35525271560688054251e-20,
twom90 = 8.07793566946316088742e-28,
twom113 = 9.62964972193617926528e-35,
twom124 = 4.70197740328915003187e-38;
-
/*
! * Extract the exponent and normalized significand (represented as
! * an array of five doubles) from a finite, nonzero quad.
! */
static int
__q_unpack(const union longdouble *x, double *s)
{
union {
double d;
unsigned int l[2];
} u;
double b;
unsigned int lx, w[3];
int ex;
/* get the normalized significand and exponent */
! ex = (int) ((x->l.msw & 0x7fffffff) >> 16);
lx = x->l.msw & 0xffff;
! if (ex)
! {
lx |= 0x10000;
w[0] = x->l.frac2;
w[1] = x->l.frac3;
w[2] = x->l.frac4;
! }
! else
! {
! if (lx | (x->l.frac2 & 0xfffe0000))
! {
w[0] = x->l.frac2;
w[1] = x->l.frac3;
w[2] = x->l.frac4;
ex = 1;
! }
! else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000))
! {
lx = x->l.frac2;
w[0] = x->l.frac3;
w[1] = x->l.frac4;
w[2] = 0;
ex = -31;
! }
! else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000))
! {
lx = x->l.frac3;
w[0] = x->l.frac4;
w[1] = w[2] = 0;
ex = -63;
! }
! else
! {
lx = x->l.frac4;
w[0] = w[1] = w[2] = 0;
ex = -95;
}
! while ((lx & 0x10000) == 0)
! {
lx = (lx << 1) | (w[0] >> 31);
w[0] = (w[0] << 1) | (w[1] >> 31);
w[1] = (w[1] << 1) | (w[2] >> 31);
w[2] <<= 1;
ex--;
--- 115,173 ----
twom66 = 1.35525271560688054251e-20,
twom90 = 8.07793566946316088742e-28,
twom113 = 9.62964972193617926528e-35,
twom124 = 4.70197740328915003187e-38;
/*
! * Extract the exponent and normalized significand (represented as
! * an array of five doubles) from a finite, nonzero quad.
! */
static int
__q_unpack(const union longdouble *x, double *s)
{
union {
double d;
unsigned int l[2];
} u;
+
double b;
unsigned int lx, w[3];
int ex;
/* get the normalized significand and exponent */
! ex = (int)((x->l.msw & 0x7fffffff) >> 16);
lx = x->l.msw & 0xffff;
!
! if (ex) {
lx |= 0x10000;
w[0] = x->l.frac2;
w[1] = x->l.frac3;
w[2] = x->l.frac4;
! } else {
! if (lx | (x->l.frac2 & 0xfffe0000)) {
w[0] = x->l.frac2;
w[1] = x->l.frac3;
w[2] = x->l.frac4;
ex = 1;
! } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) {
lx = x->l.frac2;
w[0] = x->l.frac3;
w[1] = x->l.frac4;
w[2] = 0;
ex = -31;
! } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) {
lx = x->l.frac3;
w[0] = x->l.frac4;
w[1] = w[2] = 0;
ex = -63;
! } else {
lx = x->l.frac4;
w[0] = w[1] = w[2] = 0;
ex = -95;
}
!
! while ((lx & 0x10000) == 0) {
lx = (lx << 1) | (w[0] >> 31);
w[0] = (w[0] << 1) | (w[1] >> 31);
w[1] = (w[1] << 1) | (w[2] >> 31);
w[2] <<= 1;
ex--;
*** 217,243 ****
u.l[LOWORD] = 0;
b = u.d;
u.l[LOWORD] = w[2] & 0xffffff;
s[4] = u.d - b;
! return ex - 0x3fff;
}
-
/*
! * Pack an exponent and array of three doubles representing a finite,
! * nonzero number into a quad. Assume the sign is already there and
! * the rounding mode has been fudged accordingly.
! */
static void
__q_pack(const double *z, int exp, enum fp_direction_type rm,
union longdouble *x, int *inexact)
{
union {
double d;
unsigned int l[2];
} u;
double s[3], t, t2;
unsigned int msw, frac2, frac3, frac4;
/* bias exponent and strip off integer bit */
exp += 0x3fff;
--- 205,231 ----
u.l[LOWORD] = 0;
b = u.d;
u.l[LOWORD] = w[2] & 0xffffff;
s[4] = u.d - b;
! return (ex - 0x3fff);
}
/*
! * Pack an exponent and array of three doubles representing a finite,
! * nonzero number into a quad. Assume the sign is already there and
! * the rounding mode has been fudged accordingly.
! */
static void
__q_pack(const double *z, int exp, enum fp_direction_type rm,
union longdouble *x, int *inexact)
{
union {
double d;
unsigned int l[2];
} u;
+
double s[3], t, t2;
unsigned int msw, frac2, frac3, frac4;
/* bias exponent and strip off integer bit */
exp += 0x3fff;
*** 280,320 ****
* t2 will be non-negative due to rounding mode
*/
t = s[0] + s[1];
t2 = (s[0] - t) + s[1];
! if (t != zero)
! {
*inexact = 1;
/* decide whether to round the fraction up */
if (rm == fp_positive || (rm == fp_nearest && (t > twom113 ||
! (t == twom113 && (t2 != zero || frac4 & 1)))))
! {
/* round up and renormalize if necessary */
! if (++frac4 == 0)
! if (++frac3 == 0)
! if (++frac2 == 0)
! if (++msw == 0x10000)
! {
msw = 0;
exp++;
}
}
}
/* assemble the result */
x->l.msw |= msw | (exp << 16);
x->l.frac2 = frac2;
x->l.frac3 = frac3;
x->l.frac4 = frac4;
}
-
/*
! * Compute the square root of x and place the TP result in s.
! */
static void
__q_tp_sqrt(const double *x, double *s)
{
double c, rr, r[3], tt[3], t[5];
--- 268,307 ----
* t2 will be non-negative due to rounding mode
*/
t = s[0] + s[1];
t2 = (s[0] - t) + s[1];
! if (t != zero) {
*inexact = 1;
/* decide whether to round the fraction up */
if (rm == fp_positive || (rm == fp_nearest && (t > twom113 ||
! (t == twom113 && (t2 != zero || frac4 & 1))))) {
/* round up and renormalize if necessary */
! if (++frac4 == 0) {
! if (++frac3 == 0) {
! if (++frac2 == 0) {
! if (++msw == 0x10000) {
msw = 0;
exp++;
}
}
}
+ }
+ }
+ }
/* assemble the result */
x->l.msw |= msw | (exp << 16);
x->l.frac2 = frac2;
x->l.frac3 = frac3;
x->l.frac4 = frac4;
}
/*
! * Compute the square root of x and place the TP result in s.
! */
static void
__q_tp_sqrt(const double *x, double *s)
{
double c, rr, r[3], tt[3], t[5];
*** 351,371 ****
r[1] = (r[1] - c) - t[3] * t[3];
t[4] = (rr * (r[0] + r[1]) + twom66) - twom66;
/* here we just need to get the sign of the remainder */
! c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1])
! - tt[2] * t[4]) - (t[3] + t[3]) * t[4]) - t[4] * t[4];
/* reduce to three doubles */
t[0] += t[1];
t[1] = t[2] + t[3];
t[2] = t[4];
/* if the third term might lie on a rounding boundary, perturb it */
! if (c != zero && t[2] == (twom62 + t[2]) - twom62)
! {
if (c < zero)
t[2] -= twom124;
else
t[2] += twom124;
}
--- 338,357 ----
r[1] = (r[1] - c) - t[3] * t[3];
t[4] = (rr * (r[0] + r[1]) + twom66) - twom66;
/* here we just need to get the sign of the remainder */
! c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1]) - tt[2] * t[4]) -
! (t[3] + t[3]) * t[4]) - t[4] * t[4];
/* reduce to three doubles */
t[0] += t[1];
t[1] = t[2] + t[3];
t[2] = t[4];
/* if the third term might lie on a rounding boundary, perturb it */
! if (c != zero && t[2] == (twom62 + t[2]) - twom62) {
if (c < zero)
t[2] -= twom124;
else
t[2] += twom124;
}
*** 375,400 ****
t[2] += (t[1] - c);
t[1] = c;
c = t[0] + t[1];
s[1] = t[1] + (t[0] - c);
s[0] = c;
! if (s[1] == zero)
! {
c = s[0] + t[2];
s[1] = t[2] + (s[0] - c);
s[0] = c;
s[2] = zero;
! }
! else
! {
c = s[1] + t[2];
s[2] = t[2] + (s[1] - c);
s[1] = c;
}
}
-
long double
sqrtl(long double ldx)
{
union longdouble x;
volatile double t;
--- 361,383 ----
t[2] += (t[1] - c);
t[1] = c;
c = t[0] + t[1];
s[1] = t[1] + (t[0] - c);
s[0] = c;
!
! if (s[1] == zero) {
c = s[0] + t[2];
s[1] = t[2] + (s[0] - c);
s[0] = c;
s[2] = zero;
! } else {
c = s[1] + t[2];
s[2] = t[2] + (s[1] - c);
s[1] = c;
}
}
long double
sqrtl(long double ldx)
{
union longdouble x;
volatile double t;
*** 406,479 ****
t = zero;
t -= zero;
/* check for zero operand */
x.d = ldx;
if (!((x.l.msw & 0x7fffffff) | x.l.frac2 | x.l.frac3 | x.l.frac4))
! return ldx;
/* handle nan and inf cases */
! if ((x.l.msw & 0x7fffffff) >= 0x7fff0000)
! {
! if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4)
! {
! if (!(x.l.msw & 0x8000))
! {
/* snan, signal invalid */
t += snan.d;
}
x.l.msw |= 0x8000;
! return x.d;
}
! if (x.l.msw & 0x80000000)
! {
/* sqrt(-inf), signal invalid */
t = -one;
t = sqrt(t);
! return qnan.d;
}
/* sqrt(inf), return inf */
! return x.d;
}
/* handle negative numbers */
! if (x.l.msw & 0x80000000)
! {
t = -one;
t = sqrt(t);
! return qnan.d;
}
/* now x is finite, positive */
traps = __swapTE(0);
exc = __swapEX(0);
rm = __swapRD(fp_nearest);
ex = __q_unpack(&x, xx);
! if (ex & 1)
! {
/* make exponent even */
xx[0] += xx[0];
xx[1] += xx[1];
xx[2] += xx[2];
xx[3] += xx[3];
xx[4] += xx[4];
ex--;
}
__q_tp_sqrt(xx, zz);
/* put everything together */
x.l.msw = 0;
inexact = 0;
__q_pack(zz, ex >> 1, rm, &x, &inexact);
(void) __swapRD(rm);
(void) __swapEX(exc);
(void) __swapTE(traps);
! if (inexact)
! {
t = huge;
t += tiny;
}
! return x.d;
}
--- 389,463 ----
t = zero;
t -= zero;
/* check for zero operand */
x.d = ldx;
+
if (!((x.l.msw & 0x7fffffff) | x.l.frac2 | x.l.frac3 | x.l.frac4))
! return (ldx);
/* handle nan and inf cases */
! if ((x.l.msw & 0x7fffffff) >= 0x7fff0000) {
! if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4) {
! if (!(x.l.msw & 0x8000)) {
/* snan, signal invalid */
t += snan.d;
}
+
x.l.msw |= 0x8000;
! return (x.d);
}
!
! if (x.l.msw & 0x80000000) {
/* sqrt(-inf), signal invalid */
t = -one;
t = sqrt(t);
! return (qnan.d);
}
+
/* sqrt(inf), return inf */
! return (x.d);
}
/* handle negative numbers */
! if (x.l.msw & 0x80000000) {
t = -one;
t = sqrt(t);
! return (qnan.d);
}
/* now x is finite, positive */
traps = __swapTE(0);
exc = __swapEX(0);
rm = __swapRD(fp_nearest);
ex = __q_unpack(&x, xx);
!
! if (ex & 1) {
/* make exponent even */
xx[0] += xx[0];
xx[1] += xx[1];
xx[2] += xx[2];
xx[3] += xx[3];
xx[4] += xx[4];
ex--;
}
+
__q_tp_sqrt(xx, zz);
/* put everything together */
x.l.msw = 0;
inexact = 0;
__q_pack(zz, ex >> 1, rm, &x, &inexact);
(void) __swapRD(rm);
(void) __swapEX(exc);
(void) __swapTE(traps);
!
! if (inexact) {
t = huge;
t += tiny;
}
!
! return (x.d);
}