From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 146 ------------------------------------ libs/musl/src/math/log10.c | 8 +- libs/musl/src/math/log10f.c | 6 +- 3 files changed, 9 insertions(+), 151 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index c35d997ef52..425eccb9c9e 100644 --- a/dlls/msvcrt/math.c +++ b/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. */ diff --git a/libs/musl/src/math/log10.c b/libs/musl/src/math/log10.c index df7d591c292..6ec431bfe22 100644 --- a/libs/musl/src/math/log10.c +++ b/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; diff --git a/libs/musl/src/math/log10f.c b/libs/musl/src/math/log10f.c index 64cf8c85be0..0233e2ff332 100644 --- a/libs/musl/src/math/log10f.c +++ b/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;