ntdll: Copy sqrt() implementation from msvcrt.

dlls/ntdll/math.c

@@ -1047,6 +1047,37 @@ static inline int pow_checkint(UINT64 iy)
     return 2;
 }
 
+/* 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.@)
@@ -1756,10 +1787,109 @@ double CDECL sin( double x )
 
 /*********************************************************************
  *                  sqrt   (NTDLL.@)
+ *
+ * Copied from musl: src/math/sqrt.c
  */
-double CDECL sqrt( double d )
+double CDECL sqrt( double x )
 {
-    return unix_funcs->sqrt( d );
+    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;
 }
 
 /*********************************************************************