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 ----