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); }