diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 60d85bc7be7..0f9a415d88c 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -5496,150 +5496,6 @@ void __cdecl __libm_sse2_sqrt_precise(void) #if _MSVCR_VER>=120 -/********************************************************************* - * log1p (MSVCR120.@) - * - * Copied from musl: src/math/log1p.c - */ -double CDECL log1p(double x) -{ - static const double ln2_hi = 6.93147180369123816490e-01, - ln2_lo = 1.90821492927058770002e-10, - 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, c, s, z, R, w, t1, t2, dk; - UINT32 hx, hu; - int k; - - hx = u.i >> 32; - k = 1; - if (hx < 0x3fda827a || hx >> 31) { /* 1+x < sqrt(2)+ */ - if (hx >= 0xbff00000) { /* x <= -1.0 */ - if (x == -1) { - *_errno() = ERANGE; - return x / 0.0; /* og1p(-1) = -inf */ - } - *_errno() = EDOM; - return (x-x) / 0.0; /* log1p(x<-1) = NaN */ - } - if (hx << 1 < 0x3ca00000 << 1) { /* |x| < 2**-53 */ - fp_barrier(x + 0x1p120f); - /* underflow if subnormal */ - if ((hx & 0x7ff00000) == 0) - fp_barrierf(x); - return x; - } - if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ - k = 0; - c = 0; - f = x; - } - } else if (hx >= 0x7ff00000) - return x; - if (k) { - u.f = 1 + x; - hu = u.i >> 32; - hu += 0x3ff00000 - 0x3fe6a09e; - k = (int)(hu >> 20) - 0x3ff; - /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ - if (k < 54) { - c = k >= 2 ? 1 - (u.f - x) : x - (u.f - 1); - c /= u.f; - } else - c = 0; - /* reduce u into [sqrt(2)/2, sqrt(2)] */ - hu = (hu & 0x000fffff) + 0x3fe6a09e; - u.i = (UINT64)hu << 32 | (u.i & 0xffffffff); - f = u.f - 1; - } - 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; - dk = k; - return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi; -} - -/********************************************************************* - * log1pf (MSVCR120.@) - * - * Copied from musl: src/math/log1pf.c - */ -float CDECL log1pf(float x) -{ - static const float ln2_hi = 6.9313812256e-01, - ln2_lo = 9.0580006145e-06, - 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, c, s, z, R, w, t1, t2, dk; - UINT32 ix, iu; - int k; - - ix = u.i; - k = 1; - if (ix < 0x3ed413d0 || ix >> 31) { /* 1+x < sqrt(2)+ */ - if (ix >= 0xbf800000) { /* x <= -1.0 */ - if (x == -1) { - *_errno() = ERANGE; - return x / 0.0f; /* log1p(-1)=+inf */ - } - *_errno() = EDOM; - return (x - x) / 0.0f; /* log1p(x<-1)=NaN */ - } - if (ix<<1 < 0x33800000<<1) { /* |x| < 2**-24 */ - /* underflow if subnormal */ - if ((ix & 0x7f800000) == 0) - fp_barrierf(x * x); - return x; - } - if (ix <= 0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ - k = 0; - c = 0; - f = x; - } - } else if (ix >= 0x7f800000) - return x; - if (k) { - u.f = 1 + x; - iu = u.i; - iu += 0x3f800000 - 0x3f3504f3; - k = (int)(iu >> 23) - 0x7f; - /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ - if (k < 25) { - c = k >= 2 ? 1 - (u.f - x) : x - (u.f - 1); - c /= u.f; - } else - c = 0; - /* reduce u into [sqrt(2)/2, sqrt(2)] */ - iu = (iu & 0x007fffff) + 0x3f3504f3; - u.i = iu; - f = u.f - 1; - } - 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; - dk = k; - return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi; -} - /********************************************************************* * log2 (MSVCR120.@) * diff --git a/libs/musl/src/math/log1p.c b/libs/musl/src/math/log1p.c index 9376b13ed6a..e206f17ad5e 100644 --- a/libs/musl/src/math/log1p.c +++ b/libs/musl/src/math/log1p.c @@ -77,11 +77,15 @@ double __cdecl log1p(double x) k = 1; if (hx < 0x3fda827a || hx>>31) { /* 1+x < sqrt(2)+ */ if (hx >= 0xbff00000) { /* x <= -1.0 */ - if (x == -1) + if (x == -1) { + errno = ERANGE; return x/0.0; /* log1p(-1) = -inf */ + } + errno = EDOM; return (x-x)/0.0; /* log1p(x<-1) = NaN */ } if (hx<<1 < 0x3ca00000<<1) { /* |x| < 2**-53 */ + fp_barrier(x + 0x1p120f); /* underflow if subnormal */ if ((hx&0x7ff00000) == 0) FORCE_EVAL((float)x); diff --git a/libs/musl/src/math/log1pf.c b/libs/musl/src/math/log1pf.c index e2cda58438b..9637f672e20 100644 --- a/libs/musl/src/math/log1pf.c +++ b/libs/musl/src/math/log1pf.c @@ -32,8 +32,11 @@ float __cdecl log1pf(float x) k = 1; if (ix < 0x3ed413d0 || ix>>31) { /* 1+x < sqrt(2)+ */ if (ix >= 0xbf800000) { /* x <= -1.0 */ - if (x == -1) + if (x == -1) { + errno = ERANGE; return x/0.0f; /* log1p(-1)=+inf */ + } + errno = EDOM; return (x-x)/0.0f; /* log1p(x<-1)=NaN */ } if (ix<<1 < 0x33800000<<1) { /* |x| < 2**-24 */