ntdll: Use the sqrt() implementation from the bundled musl library.

dlls/ntdll/math.c

@@ -47,38 +47,6 @@ double math_error( int type, const char *name, double arg1, double arg2, double
     return retval;
 }
 
-/* Copied from musl: src/math/__fpclassify.c */
-static short _dclass(double x)
-{
-    union { double f; UINT64 i; } u = { x };
-    int e = u.i >> 52 & 0x7ff;
-
-    if (!e) return u.i << 1 ? FP_SUBNORMAL : FP_ZERO;
-    if (e == 0x7ff) return (u.i << 12) ? FP_NAN : FP_INFINITE;
-    return FP_NORMAL;
-}
-
-static BOOL sqrt_validate( double *x, BOOL update_sw )
-{
-    short c = _dclass(*x);
-
-    if (c == FP_ZERO) return FALSE;
-    if (c == FP_NAN)
-    {
-        /* set signaling bit */
-        *(ULONGLONG*)x |= 0x8000000000000ULL;
-        return FALSE;
-    }
-    if (signbit(*x))
-    {
-        *x = -NAN;
-        return FALSE;
-    }
-    if (c == FP_INFINITE) return FALSE;
-    return TRUE;
-}
-
-
 /*********************************************************************
  *                  abs   (NTDLL.@)
  */
@@ -157,113 +125,6 @@ double CDECL floor( double x )
     return u.f;
 }
 
-/*********************************************************************
- *                  sqrt   (NTDLL.@)
- *
- * Copied from musl: src/math/sqrt.c
- */
-double CDECL sqrt( double x )
-{
-    static const double tiny = 1.0e-300;
-
-    double z;
-    int sign = 0x80000000;
-    int ix0,s0,q,m,t,i;
-    unsigned int r,t1,s1,ix1,q1;
-    ULONGLONG ix;
-
-    if (!sqrt_validate(&x, TRUE))
-        return x;
-
-    ix = *(ULONGLONG*)&x;
-    ix0 = ix >> 32;
-    ix1 = ix;
-
-    /* normalize x */
-    m = ix0 >> 20;
-    if (m == 0) {  /* subnormal x */
-        while (ix0 == 0) {
-            m -= 21;
-            ix0 |= (ix1 >> 11);
-            ix1 <<= 21;
-        }
-        for (i=0; (ix0 & 0x00100000) == 0; i++)
-            ix0 <<= 1;
-        m -= i - 1;
-        ix0 |= ix1 >> (32 - i);
-        ix1 <<= i;
-    }
-    m -= 1023;    /* unbias exponent */
-    ix0 = (ix0 & 0x000fffff) | 0x00100000;
-    if (m & 1) {  /* odd m, double x to make it even */
-        ix0 += ix0 + ((ix1 & sign) >> 31);
-        ix1 += ix1;
-    }
-    m >>= 1;      /* m = [m/2] */
-
-    /* generate sqrt(x) bit by bit */
-    ix0 += ix0 + ((ix1 & sign) >> 31);
-    ix1 += ix1;
-    q = q1 = s0 = s1 = 0;  /* [q,q1] = sqrt(x) */
-    r = 0x00200000;        /* r = moving bit from right to left */
-
-    while (r != 0) {
-        t = s0 + r;
-        if (t <= ix0) {
-            s0   = t + r;
-            ix0 -= t;
-            q   += r;
-        }
-        ix0 += ix0 + ((ix1 & sign) >> 31);
-        ix1 += ix1;
-        r >>= 1;
-    }
-
-    r = sign;
-    while (r != 0) {
-        t1 = s1 + r;
-        t  = s0;
-        if (t < ix0 || (t == ix0 && t1 <= ix1)) {
-            s1 = t1 + r;
-            if ((t1&sign) == sign && (s1 & sign) == 0)
-                s0++;
-            ix0 -= t;
-            if (ix1 < t1)
-                ix0--;
-            ix1 -= t1;
-            q1 += r;
-        }
-        ix0 += ix0 + ((ix1 & sign) >> 31);
-        ix1 += ix1;
-        r >>= 1;
-    }
-
-    /* use floating add to find out rounding direction */
-    if ((ix0 | ix1) != 0) {
-        z = 1.0 - tiny; /* raise inexact flag */
-        if (z >= 1.0) {
-            z = 1.0 + tiny;
-            if (q1 == (unsigned int)0xffffffff) {
-                q1 = 0;
-                q++;
-            } else if (z > 1.0) {
-                if (q1 == (unsigned int)0xfffffffe)
-                    q++;
-                q1 += 2;
-            } else
-                q1 += q1 & 1;
-        }
-    }
-    ix0 = (q >> 1) + 0x3fe00000;
-    ix1 = q1 >> 1;
-    if (q & 1)
-        ix1 |= sign;
-    ix = ix0 + ((unsigned int)m << 20);
-    ix <<= 32;
-    ix |= ix1;
-    return *(double*)&ix;
-}
-
 #if (defined(__GNUC__) || defined(__clang__)) && defined(__i386__)
 
 #define FPU_DOUBLE(var) double var; \