Print this page
5262 libm needs to be carefully unifdef'd
5268 libm doesn't need to hide symbols which are already local

*** 28,38 **** #pragma weak sqrt = __sqrt #include "libm.h" - #ifdef __INLINE extern double __inline_sqrt(double); double sqrt(double x) { --- 28,37 ----
*** 41,150 **** if (isnan(x)) return (z); return ((x < 0.0)? _SVID_libm_err(x, x, 26) : z); } - #else /* defined(__INLINE) */ - - /* - * Warning: This correctly rounded sqrt is extremely slow because it computes - * the sqrt bit by bit using integer arithmetic. - */ - - static const double big = 1.0e30, small = 1.0e-30; - - double - sqrt(double x) - { - double z; - unsigned r, t1, s1, ix1, q1; - int ix0, s0, j, q, m, n, t; - int *px = (int *)&x, *pz = (int *)&z; - - ix0 = px[HIWORD]; - ix1 = px[LOWORD]; - if ((ix0 & 0x7ff00000) == 0x7ff00000) { /* x is inf or NaN */ - if (ix0 == 0xfff00000 && ix1 == 0) - return (_SVID_libm_err(x, x, 26)); - return (x + x); - } - if (((ix0 & 0x7fffffff) | ix1) == 0) /* x is zero */ - return (x); - - /* extract exponent and significand */ - m = ilogb(x); - z = scalbn(x, -m); - ix0 = (pz[HIWORD] & 0x000fffff) | 0x00100000; - ix1 = pz[LOWORD]; - n = m >> 1; - if (n + n != m) { - ix0 = (ix0 << 1) | (ix1 >> 31); - ix1 <<= 1; - m -= 1; - } - - /* generate sqrt(x) bit by bit */ - ix0 = (ix0 << 1) | (ix1 >> 31); - ix1 <<= 1; - q = q1 = s0 = s1 = 0; - r = 0x00200000; - - for (j = 1; j <= 22; j++) { - t = s0 + r; - if (t <= ix0) { - s0 = t + r; - ix0 -= t; - q += r; - } - ix0 = (ix0 << 1) | (ix1 >> 31); - ix1 <<= 1; - r >>= 1; - } - - r = 0x80000000; - for (j = 1; j <= 32; j++) { - t1 = s1 + r; - t = s0; - if (t < ix0 || (t == ix0 && t1 <= ix1)) { - s1 = t1 + r; - if ((t1 & 0x80000000) == 0x80000000 && - (s1 & 0x80000000) == 0) - s0 += 1; - ix0 -= t; - if (ix1 < t1) - ix0 -= 1; - ix1 -= t1; - q1 += r; - } - ix0 = (ix0 << 1) | (ix1 >> 31); - ix1 <<= 1; - r >>= 1; - } - - /* round */ - if ((ix0 | ix1) == 0) - goto done; - z = big - small; /* trigger inexact flag */ - if (z < big) - goto done; - if (q1 == 0xffffffff) { - q1 = 0; - q += 1; - goto done; - } - z = big + small; - if (z > big) { - if (q1 == 0xfffffffe) - q += 1; - q1 += 2; - goto done; - } - q1 += (q1 & 1); - done: - pz[HIWORD] = (q >> 1) + 0x3fe00000; - pz[LOWORD] = q1 >> 1; - if ((q & 1) == 1) - pz[LOWORD] |= 0x80000000; - return (scalbn(z, n)); - } - - #endif /* defined(__INLINE) */ --- 40,44 ----