msvcrt: Use the log10()/log10f() implementation from the bundled musl library.

dlls/msvcrt/math.c

@@ -1084,69 +1084,6 @@ float CDECL expf( float x )
     return y;
 }
 
-/*********************************************************************
- *      log10f (MSVCRT.@)
- */
-float CDECL log10f( float x )
-{
-    static const float ivln10hi = 4.3432617188e-01,
-        ivln10lo = -3.1689971365e-05,
-        log10_2hi = 3.0102920532e-01,
-        log10_2lo = 7.9034151668e-07,
-        Lg1 = 0xaaaaaa.0p-24,
-        Lg2 = 0xccce13.0p-25,
-        Lg3 = 0x91e9ee.0p-25,
-        Lg4 = 0xf89e26.0p-26;
-
-    union {float f; UINT32 i;} u = {x};
-    float hfsq, f, s, z, R, w, t1, t2, dk, hi, lo;
-    UINT32 ix;
-    int k;
-
-    ix = u.i;
-    k = 0;
-    if (ix < 0x00800000 || ix >> 31) { /* x < 2**-126 */
-        if (ix << 1 == 0)
-            return math_error(_SING, "log10f", x, 0, -1 / (x * x));
-        if ((ix & ~(1u << 31)) > 0x7f800000)
-            return x;
-        if (ix >> 31)
-            return math_error(_DOMAIN, "log10f", x, 0, (x - x) / (x - x));
-        /* subnormal number, scale up x */
-        k -= 25;
-        x *= 0x1p25f;
-        u.f = x;
-        ix = u.i;
-    } else if (ix >= 0x7f800000) {
-        return x;
-    } else if (ix == 0x3f800000)
-        return 0;
-
-    /* reduce x into [sqrt(2)/2, sqrt(2)] */
-    ix += 0x3f800000 - 0x3f3504f3;
-    k += (int)(ix >> 23) - 0x7f;
-    ix = (ix & 0x007fffff) + 0x3f3504f3;
-    u.i = ix;
-    x = u.f;
-
-    f = x - 1.0f;
-    s = f / (2.0f + f);
-    z = s * s;
-    w = z * z;
-    t1= w * (Lg2 + w * Lg4);
-    t2= z * (Lg1 + w * Lg3);
-    R = t2 + t1;
-    hfsq = 0.5f * f * f;
-
-    hi = f - hfsq;
-    u.f = hi;
-    u.i &= 0xfffff000;
-    hi = u.f;
-    lo = f - hi - hfsq + s * (hfsq + R);
-    dk = k;
-    return dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi + hi * ivln10hi + dk * log10_2hi;
-}
-
 /* Subnormal input is normalized so ix has negative biased exponent.
    Output is multiplied by POWF_SCALE (where 1 << 5). */
 static double powf_log2(UINT32 ix)
@@ -2599,89 +2536,6 @@ double CDECL exp( double x )
     return scale + scale * tmp;
 }
 
-/*********************************************************************
- *		log10 (MSVCRT.@)
- */
-double CDECL log10( double x )
-{
-    static const double ivln10hi = 4.34294481878168880939e-01,
-        ivln10lo = 2.50829467116452752298e-11,
-        log10_2hi = 3.01029995663611771306e-01,
-        log10_2lo = 3.69423907715893078616e-13,
-        Lg1 = 6.666666666666735130e-01,
-        Lg2 = 3.999999999940941908e-01,
-        Lg3 = 2.857142874366239149e-01,
-        Lg4 = 2.222219843214978396e-01,
-        Lg5 = 1.818357216161805012e-01,
-        Lg6 = 1.531383769920937332e-01,
-        Lg7 = 1.479819860511658591e-01;
-
-    union {double f; UINT64 i;} u = {x};
-    double hfsq, f, s, z, R, w, t1, t2, dk, y, hi, lo, val_hi, val_lo;
-    UINT32 hx;
-    int k;
-
-    hx = u.i >> 32;
-    k = 0;
-    if (hx < 0x00100000 || hx >> 31) {
-        if (u.i << 1 == 0)
-            return math_error(_SING, "log10", x, 0, -1 / (x * x));
-        if ((u.i & ~(1ULL << 63)) > 0x7ff0000000000000ULL)
-            return x;
-        if (hx >> 31)
-            return math_error(_DOMAIN, "log10", x, 0, (x - x) / (x - x));
-        /* subnormal number, scale x up */
-        k -= 54;
-        x *= 0x1p54;
-        u.f = x;
-        hx = u.i >> 32;
-    } else if (hx >= 0x7ff00000) {
-        return x;
-    } else if (hx == 0x3ff00000 && u.i<<32 == 0)
-        return 0;
-
-    /* reduce x into [sqrt(2)/2, sqrt(2)] */
-    hx += 0x3ff00000 - 0x3fe6a09e;
-    k += (int)(hx >> 20) - 0x3ff;
-    hx = (hx & 0x000fffff) + 0x3fe6a09e;
-    u.i = (UINT64)hx << 32 | (u.i & 0xffffffff);
-    x = u.f;
-
-    f = x - 1.0;
-    hfsq = 0.5 * f * f;
-    s = f / (2.0 + f);
-    z = s * s;
-    w = z * z;
-    t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
-    t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
-    R = t2 + t1;
-
-    /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
-    hi = f - hfsq;
-    u.f = hi;
-    u.i &= (UINT64)-1 << 32;
-    hi = u.f;
-    lo = f - hi - hfsq + s * (hfsq + R);
-
-    /* val_hi+val_lo ~ log10(1+f) + k*log10(2) */
-    val_hi = hi * ivln10hi;
-    dk = k;
-    y = dk * log10_2hi;
-    val_lo = dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi;
-
-    /*
-     * Extra precision in for adding y is not strictly needed
-     * since there is no very large cancellation near x = sqrt(2) or
-     * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
-     * with some parallelism and it reduces the error for many args.
-     */
-    w = y + val_hi;
-    val_lo += (y - w) + val_hi;
-    val_hi = w;
-
-    return val_lo + val_hi;
-}
-
 /* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
    additional 15 bits precision. IX is the bit representation of x, but
    normalized in the subnormal range using the sign bit for the exponent. */

libs/musl/src/math/log10.c

@@ -45,9 +45,11 @@ double __cdecl log10(double x)
 	k = 0;
 	if (hx < 0x00100000 || hx>>31) {
 		if (u.i<<1 == 0)
-			return -1/(x*x);  /* log(+-0)=-inf */
-		if (hx>>31)
-			return (x-x)/0.0; /* log(-#) = NaN */
+			return math_error(_SING, "log10", x, 0, -1 / (x * x));
+		if ((u.i & ~(1ULL << 63)) > 0x7ff0000000000000ULL)
+			return x;
+		if (hx >> 31)
+			return math_error(_DOMAIN, "log10", x, 0, (x - x) / (x - x));
 		/* subnormal number, scale x up */
 		k -= 54;
 		x *= 0x1p54;

libs/musl/src/math/log10f.c

@@ -39,9 +39,11 @@ float __cdecl log10f(float x)
 	k = 0;
 	if (ix < 0x00800000 || ix>>31) {  /* x < 2**-126  */
 		if (ix<<1 == 0)
-			return -1/(x*x);  /* log(+-0)=-inf */
+			return math_error(_SING, "log10f", x, 0, -1 / (x * x));
+		if ((ix & ~(1u << 31)) > 0x7f800000)
+			return x;
 		if (ix>>31)
-			return (x-x)/0.0f; /* log(-#) = NaN */
+			return math_error(_DOMAIN, "log10f", x, 0, (x - x) / (x - x));
 		/* subnormal number, scale up x */
 		k -= 25;
 		x *= 0x1p25f;