fix double-processing of DT_RELR relocations in ldso relocating itself

[musl] / src / math / cbrtl.c

diff --git a/src/math/cbrtl.c b/src/math/cbrtl.c
index 5297d68..ceff913 100644 (file)
--- a/src/math/cbrtl.c
+++ b/src/math/cbrtl.c
@@ -23,73 +23,50 @@ long double cbrtl(long double x)
        return cbrt(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#define BIAS    (LDBL_MAX_EXP - 1)
-static const unsigned
-B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
+static const unsigned B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
 
 long double cbrtl(long double x)
 {
-       union IEEEl2bits u, v;
+       union ldshape u = {x}, v;
+       union {float f; uint32_t i;} uft;
        long double r, s, t, w;
-       double dr, dt, dx;
-       float ft, fx;
-       uint32_t hx;
-       uint16_t expsign;
-       int k;
-
-       u.e = x;
-       expsign = u.xbits.expsign;
-       k = expsign & 0x7fff;
+       double_t dr, dt, dx;
+       float_t ft;
+       int e = u.i.se & 0x7fff;
+       int sign = u.i.se & 0x8000;
 
        /*
         * If x = +-Inf, then cbrt(x) = +-Inf.
         * If x = NaN, then cbrt(x) = NaN.
         */
-       if (k == BIAS + LDBL_MAX_EXP)
+       if (e == 0x7fff)
                return x + x;
-
-// FIXME: extended precision is default on linux..
-#undef __i386__
-#ifdef __i386__
-       fp_prec_t oprec;
-
-       oprec = fpgetprec();
-       if (oprec != FP_PE)
-               fpsetprec(FP_PE);
-#endif
-
-       if (k == 0) {
-               /* If x = +-0, then cbrt(x) = +-0. */
-               if ((u.bits.manh | u.bits.manl) == 0) {
-#ifdef __i386__
-                       if (oprec != FP_PE)
-                               fpsetprec(oprec);
-#endif
-                       return (x);
-               }
+       if (e == 0) {
                /* Adjust subnormal numbers. */
-               u.e *= 0x1.0p514;
-               k = u.bits.exp;
-               k -= BIAS + 514;
-       } else
-               k -= BIAS;
-       u.xbits.expsign = BIAS;
-       v.e = 1;
-
-       x = u.e;
-       switch (k % 3) {
+               u.f *= 0x1p120;
+               e = u.i.se & 0x7fff;
+               /* If x = +-0, then cbrt(x) = +-0. */
+               if (e == 0)
+                       return x;
+               e -= 120;
+       }
+       e -= 0x3fff;
+       u.i.se = 0x3fff;
+       x = u.f;
+       switch (e % 3) {
        case 1:
        case -2:
-               x = 2*x;
-               k--;
+               x *= 2;
+               e--;
                break;
        case 2:
        case -1:
-               x = 4*x;
-               k -= 2;
+               x *= 4;
+               e -= 2;
                break;
        }
-       v.xbits.expsign = (expsign & 0x8000) | (BIAS + k / 3);
+       v.f = 1.0;
+       v.i.se = sign | (0x3fff + e/3);
 
        /*
         * The following is the guts of s_cbrtf, with the handling of
@@ -98,9 +75,9 @@ long double cbrtl(long double x)
         */
 
        /* ~5-bit estimate: */
-       fx = x;
-       GET_FLOAT_WORD(hx, fx);
-       SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1));
+       uft.f = x;
+       uft.i = (uft.i & 0x7fffffff)/3 + B1;
+       ft = uft.f;
 
        /* ~16-bit estimate: */
        dx = x;
@@ -118,7 +95,7 @@ long double cbrtl(long double x)
         * Round it away from zero to 32 bits (32 so that t*t is exact, and
         * away from zero for technical reasons).
         */
-       t = dt + (0x1.0p32L + 0x1.0p-32L) - 0x1.0p32;
+       t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32;
 #elif LDBL_MANT_DIG == 113
        /*
         * Round dt away from zero to 47 bits.  Since we don't trust the 47,
@@ -129,8 +106,6 @@ long double cbrtl(long double x)
         * dt is much more accurate than needed.
         */
        t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
-#else
-#error "Unsupported long double format"
 #endif
 
        /*
@@ -143,11 +118,7 @@ long double cbrtl(long double x)
        r = (r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
        t = t+t*r;       /* error <= 0.5 + 0.5/3 + epsilon */
 
-       t *= v.e;
-#ifdef __i386__
-       if (oprec != FP_PE)
-               fpsetprec(oprec);
-#endif
+       t *= v.f;
        return t;
 }
 #endif