From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 597 -------------------------------------- libs/musl/src/math/pow.c | 42 ++- libs/musl/src/math/powf.c | 21 +- 3 files changed, 25 insertions(+), 635 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 71d193290de..1e67ced4360 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -540,184 +540,6 @@ float CDECL expf( float x ) return y; }
-/* 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) -{ - static const struct { - double invc, logc; - } T[] = { - { 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 * (1 << 5) }, - { 0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2 * (1 << 5) }, - { 0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2 * (1 << 5) }, - { 0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2 * (1 << 5) }, - { 0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2 * (1 << 5) }, - { 0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3 * (1 << 5) }, - { 0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3 * (1 << 5) }, - { 0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4 * (1 << 5) }, - { 0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5 * (1 << 5) }, - { 0x1p+0, 0x0p+0 * (1 << 4) }, - { 0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4 * (1 << 5) }, - { 0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3 * (1 << 5) }, - { 0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3 * (1 << 5) }, - { 0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2 * (1 << 5) }, - { 0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2 * (1 << 5) }, - { 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 * (1 << 5) } - }; - static const double A[] = { - 0x1.27616c9496e0bp-2 * (1 << 5), -0x1.71969a075c67ap-2 * (1 << 5), - 0x1.ec70a6ca7baddp-2 * (1 << 5), -0x1.7154748bef6c8p-1 * (1 << 5), - 0x1.71547652ab82bp0 * (1 << 5) - }; - - double z, r, r2, r4, p, q, y, y0, invc, logc; - UINT32 iz, top, tmp; - int k, i; - - /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. - The range is split into N subintervals. - The ith subinterval contains z and c is near its center. */ - tmp = ix - 0x3f330000; - i = (tmp >> (23 - 4)) % (1 << 4); - top = tmp & 0xff800000; - iz = ix - top; - k = (INT32)top >> (23 - 5); /* arithmetic shift */ - invc = T[i].invc; - logc = T[i].logc; - z = *(float*)&iz; - - /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ - r = z * invc - 1; - y0 = logc + (double)k; - - /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ - r2 = r * r; - y = A[0] * r + A[1]; - p = A[2] * r + A[3]; - r4 = r2 * r2; - q = A[4] * r + y0; - q = p * r2 + q; - y = y * r4 + q; - return y; -} - -/* The output of log2 and thus the input of exp2 is either scaled by N - (in case of fast toint intrinsics) or not. The unscaled xd must be - in [-1021,1023], sign_bias sets the sign of the result. */ -static float powf_exp2(double xd, UINT32 sign_bias) -{ - static const double C[] = { - 0x1.c6af84b912394p-5 / (1 << 5) / (1 << 5) / (1 << 5), - 0x1.ebfce50fac4f3p-3 / (1 << 5) / (1 << 5), - 0x1.62e42ff0c52d6p-1 / (1 << 5) - }; - - UINT64 ki, ski, t; - double kd, z, r, r2, y, s; - - /* N*x = k + r with r in [-1/2, 1/2] */ - kd = round(xd); /* k */ - ki = (INT64)kd; - r = xd - kd; - - /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ - t = exp2f_T[ki % (1 << 5)]; - ski = ki + sign_bias; - t += ski << (52 - 5); - s = *(double*)&t; - z = C[0] * r + C[1]; - r2 = r * r; - y = C[2] * r + 1; - y = z * r2 + y; - y = y * s; - return y; -} - -/* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is - the bit representation of a non-zero finite floating-point value. */ -static int powf_checkint(UINT32 iy) -{ - int e = iy >> 23 & 0xff; - if (e < 0x7f) - return 0; - if (e > 0x7f + 23) - return 2; - if (iy & ((1 << (0x7f + 23 - e)) - 1)) - return 0; - if (iy & (1 << (0x7f + 23 - e))) - return 1; - return 2; -} - -/********************************************************************* - * powf (MSVCRT.@) - * - * Copied from musl: src/math/powf.c src/math/powf_data.c - */ -float CDECL powf( float x, float y ) -{ - UINT32 sign_bias = 0; - UINT32 ix, iy; - double logx, ylogx; - - ix = *(UINT32*)&x; - iy = *(UINT32*)&y; - if (ix - 0x00800000 >= 0x7f800000 - 0x00800000 || - 2 * iy - 1 >= 2u * 0x7f800000 - 1) { - /* Either (x < 0x1p-126 or inf or nan) or (y is 0 or inf or nan). */ - if (2 * iy - 1 >= 2u * 0x7f800000 - 1) { - if (2 * iy == 0) - return 1.0f; - if (ix == 0x3f800000) - return 1.0f; - if (2 * ix > 2u * 0x7f800000 || 2 * iy > 2u * 0x7f800000) - return x + y; - if (2 * ix == 2 * 0x3f800000) - return 1.0f; - if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000)) - return 0.0f; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ - return y * y; - } - if (2 * ix - 1 >= 2u * 0x7f800000 - 1) { - float x2 = x * x; - if (ix & 0x80000000 && powf_checkint(iy) == 1) - x2 = -x2; - if (iy & 0x80000000 && x2 == 0.0) - return math_error(_SING, "powf", x, y, 1 / x2); - /* Without the barrier some versions of clang hoist the 1/x2 and - thus division by zero exception can be signaled spuriously. */ - return iy & 0x80000000 ? fp_barrierf(1 / x2) : x2; - } - /* x and y are non-zero finite. */ - if (ix & 0x80000000) { - /* Finite x < 0. */ - int yint = powf_checkint(iy); - if (yint == 0) - return math_error(_DOMAIN, "powf", x, y, 0 / (x - x)); - if (yint == 1) - sign_bias = 1 << (5 + 11); - ix &= 0x7fffffff; - } - if (ix < 0x00800000) { - /* Normalize subnormal x so exponent becomes negative. */ - x *= 0x1p23f; - ix = *(UINT32*)&x; - ix &= 0x7fffffff; - ix -= 23 << 23; - } - } - logx = powf_log2(ix); - ylogx = y * logx; /* cannot overflow, y is single prec. */ - if ((*(UINT64*)&ylogx >> 47 & 0xffff) >= 0x40af800000000000llu >> 47) { - /* |y*log(x)| >= 126. */ - if (ylogx > 0x1.fffffffd1d571p+6 * (1 << 5)) - return math_error(_OVERFLOW, "powf", x, y, (sign_bias ? -1.0 : 1.0) * 0x1p1023); - if (ylogx <= -150.0 * (1 << 5)) - return math_error(_UNDERFLOW, "powf", x, y, (sign_bias ? -1.0 : 1.0) / 0x1p1023); - } - return powf_exp2(ylogx, sign_bias); -} - static BOOL sqrtf_validate( float *x ) { short c = _fdclass(*x); @@ -1408,425 +1230,6 @@ double CDECL exp( double x ) return scale + scale * tmp; }
-/* 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. */ -static double pow_log(UINT64 ix, double *tail) -{ - static const struct { - double invc, logc, logctail; - } T[] = { - {0x1.6a00000000000p+0, -0x1.62c82f2b9c800p-2, 0x1.ab42428375680p-48}, - {0x1.6800000000000p+0, -0x1.5d1bdbf580800p-2, -0x1.ca508d8e0f720p-46}, - {0x1.6600000000000p+0, -0x1.5767717455800p-2, -0x1.362a4d5b6506dp-45}, - {0x1.6400000000000p+0, -0x1.51aad872df800p-2, -0x1.684e49eb067d5p-49}, - {0x1.6200000000000p+0, -0x1.4be5f95777800p-2, -0x1.41b6993293ee0p-47}, - {0x1.6000000000000p+0, -0x1.4618bc21c6000p-2, 0x1.3d82f484c84ccp-46}, - {0x1.5e00000000000p+0, -0x1.404308686a800p-2, 0x1.c42f3ed820b3ap-50}, - {0x1.5c00000000000p+0, -0x1.3a64c55694800p-2, 0x1.0b1c686519460p-45}, - {0x1.5a00000000000p+0, -0x1.347dd9a988000p-2, 0x1.5594dd4c58092p-45}, - {0x1.5800000000000p+0, -0x1.2e8e2bae12000p-2, 0x1.67b1e99b72bd8p-45}, - {0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46}, - {0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46}, - {0x1.5400000000000p+0, -0x1.22941fbcf7800p-2, -0x1.65a242853da76p-46}, - {0x1.5200000000000p+0, -0x1.1c898c1699800p-2, -0x1.fafbc68e75404p-46}, - {0x1.5000000000000p+0, -0x1.1675cababa800p-2, 0x1.f1fc63382a8f0p-46}, - {0x1.4e00000000000p+0, -0x1.1058bf9ae4800p-2, -0x1.6a8c4fd055a66p-45}, - {0x1.4c00000000000p+0, -0x1.0a324e2739000p-2, -0x1.c6bee7ef4030ep-47}, - {0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48}, - {0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48}, - {0x1.4800000000000p+0, -0x1.fb9186d5e4000p-3, 0x1.d572aab993c87p-47}, - {0x1.4600000000000p+0, -0x1.ef0adcbdc6000p-3, 0x1.b26b79c86af24p-45}, - {0x1.4400000000000p+0, -0x1.e27076e2af000p-3, -0x1.72f4f543fff10p-46}, - {0x1.4200000000000p+0, -0x1.d5c216b4fc000p-3, 0x1.1ba91bbca681bp-45}, - {0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45}, - {0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45}, - {0x1.3e00000000000p+0, -0x1.bc286742d9000p-3, 0x1.94eb0318bb78fp-46}, - {0x1.3c00000000000p+0, -0x1.af3c94e80c000p-3, 0x1.a4e633fcd9066p-52}, - {0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45}, - {0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45}, - {0x1.3800000000000p+0, -0x1.9525a9cf45000p-3, -0x1.ad1d904c1d4e3p-45}, - {0x1.3600000000000p+0, -0x1.87fa06520d000p-3, 0x1.bbdbf7fdbfa09p-45}, - {0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45}, - {0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45}, - {0x1.3200000000000p+0, -0x1.6d60fe719d000p-3, -0x1.0e46aa3b2e266p-46}, - {0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46}, - {0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46}, - {0x1.2e00000000000p+0, -0x1.526e5e3a1b000p-3, -0x1.0de8b90075b8fp-45}, - {0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46}, - {0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46}, - {0x1.2a00000000000p+0, -0x1.371fc201e9000p-3, 0x1.178864d27543ap-48}, - {0x1.2800000000000p+0, -0x1.29552f81ff000p-3, -0x1.48d301771c408p-45}, - {0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45}, - {0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45}, - {0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47}, - {0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47}, - {0x1.2200000000000p+0, -0x1.fec9131dbe000p-4, -0x1.575545ca333f2p-45}, - {0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45}, - {0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45}, - {0x1.1e00000000000p+0, -0x1.c5e548f5bc000p-4, -0x1.d0c57585fbe06p-46}, - {0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45}, - {0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45}, - {0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46}, - {0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46}, - {0x1.1800000000000p+0, -0x1.6f0d28ae56000p-4, -0x1.69737c93373dap-45}, - {0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46}, - {0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46}, - {0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45}, - {0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45}, - {0x1.1200000000000p+0, -0x1.16536eea38000p-4, 0x1.47c5e768fa309p-46}, - {0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45}, - {0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45}, - {0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46}, - {0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46}, - {0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45}, - {0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45}, - {0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48}, - {0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48}, - {0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45}, - {0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45}, - {0x1.0600000000000p+0, -0x1.7b91b07d58000p-6, -0x1.88d5493faa639p-45}, - {0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50}, - {0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50}, - {0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46}, - {0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46}, - {0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0}, - {0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0}, - {0x1.fc00000000000p-1, 0x1.0101575890000p-7, -0x1.0c76b999d2be8p-46}, - {0x1.f800000000000p-1, 0x1.0205658938000p-6, -0x1.3dc5b06e2f7d2p-45}, - {0x1.f400000000000p-1, 0x1.8492528c90000p-6, -0x1.aa0ba325a0c34p-45}, - {0x1.f000000000000p-1, 0x1.0415d89e74000p-5, 0x1.111c05cf1d753p-47}, - {0x1.ec00000000000p-1, 0x1.466aed42e0000p-5, -0x1.c167375bdfd28p-45}, - {0x1.e800000000000p-1, 0x1.894aa149fc000p-5, -0x1.97995d05a267dp-46}, - {0x1.e400000000000p-1, 0x1.ccb73cdddc000p-5, -0x1.a68f247d82807p-46}, - {0x1.e200000000000p-1, 0x1.eea31c006c000p-5, -0x1.e113e4fc93b7bp-47}, - {0x1.de00000000000p-1, 0x1.1973bd1466000p-4, -0x1.5325d560d9e9bp-45}, - {0x1.da00000000000p-1, 0x1.3bdf5a7d1e000p-4, 0x1.cc85ea5db4ed7p-45}, - {0x1.d600000000000p-1, 0x1.5e95a4d97a000p-4, -0x1.c69063c5d1d1ep-45}, - {0x1.d400000000000p-1, 0x1.700d30aeac000p-4, 0x1.c1e8da99ded32p-49}, - {0x1.d000000000000p-1, 0x1.9335e5d594000p-4, 0x1.3115c3abd47dap-45}, - {0x1.cc00000000000p-1, 0x1.b6ac88dad6000p-4, -0x1.390802bf768e5p-46}, - {0x1.ca00000000000p-1, 0x1.c885801bc4000p-4, 0x1.646d1c65aacd3p-45}, - {0x1.c600000000000p-1, 0x1.ec739830a2000p-4, -0x1.dc068afe645e0p-45}, - {0x1.c400000000000p-1, 0x1.fe89139dbe000p-4, -0x1.534d64fa10afdp-45}, - {0x1.c000000000000p-1, 0x1.1178e8227e000p-3, 0x1.1ef78ce2d07f2p-45}, - {0x1.be00000000000p-1, 0x1.1aa2b7e23f000p-3, 0x1.ca78e44389934p-45}, - {0x1.ba00000000000p-1, 0x1.2d1610c868000p-3, 0x1.39d6ccb81b4a1p-47}, - {0x1.b800000000000p-1, 0x1.365fcb0159000p-3, 0x1.62fa8234b7289p-51}, - {0x1.b400000000000p-1, 0x1.4913d8333b000p-3, 0x1.5837954fdb678p-45}, - {0x1.b200000000000p-1, 0x1.527e5e4a1b000p-3, 0x1.633e8e5697dc7p-45}, - {0x1.ae00000000000p-1, 0x1.6574ebe8c1000p-3, 0x1.9cf8b2c3c2e78p-46}, - {0x1.ac00000000000p-1, 0x1.6f0128b757000p-3, -0x1.5118de59c21e1p-45}, - {0x1.aa00000000000p-1, 0x1.7898d85445000p-3, -0x1.c661070914305p-46}, - {0x1.a600000000000p-1, 0x1.8beafeb390000p-3, -0x1.73d54aae92cd1p-47}, - {0x1.a400000000000p-1, 0x1.95a5adcf70000p-3, 0x1.7f22858a0ff6fp-47}, - {0x1.a000000000000p-1, 0x1.a93ed3c8ae000p-3, -0x1.8724350562169p-45}, - {0x1.9e00000000000p-1, 0x1.b31d8575bd000p-3, -0x1.c358d4eace1aap-47}, - {0x1.9c00000000000p-1, 0x1.bd087383be000p-3, -0x1.d4bc4595412b6p-45}, - {0x1.9a00000000000p-1, 0x1.c6ffbc6f01000p-3, -0x1.1ec72c5962bd2p-48}, - {0x1.9600000000000p-1, 0x1.db13db0d49000p-3, -0x1.aff2af715b035p-45}, - {0x1.9400000000000p-1, 0x1.e530effe71000p-3, 0x1.212276041f430p-51}, - {0x1.9200000000000p-1, 0x1.ef5ade4dd0000p-3, -0x1.a211565bb8e11p-51}, - {0x1.9000000000000p-1, 0x1.f991c6cb3b000p-3, 0x1.bcbecca0cdf30p-46}, - {0x1.8c00000000000p-1, 0x1.07138604d5800p-2, 0x1.89cdb16ed4e91p-48}, - {0x1.8a00000000000p-1, 0x1.0c42d67616000p-2, 0x1.7188b163ceae9p-45}, - {0x1.8800000000000p-1, 0x1.1178e8227e800p-2, -0x1.c210e63a5f01cp-45}, - {0x1.8600000000000p-1, 0x1.16b5ccbacf800p-2, 0x1.b9acdf7a51681p-45}, - {0x1.8400000000000p-1, 0x1.1bf99635a6800p-2, 0x1.ca6ed5147bdb7p-45}, - {0x1.8200000000000p-1, 0x1.214456d0eb800p-2, 0x1.a87deba46baeap-47}, - {0x1.7e00000000000p-1, 0x1.2bef07cdc9000p-2, 0x1.a9cfa4a5004f4p-45}, - {0x1.7c00000000000p-1, 0x1.314f1e1d36000p-2, -0x1.8e27ad3213cb8p-45}, - {0x1.7a00000000000p-1, 0x1.36b6776be1000p-2, 0x1.16ecdb0f177c8p-46}, - {0x1.7800000000000p-1, 0x1.3c25277333000p-2, 0x1.83b54b606bd5cp-46}, - {0x1.7600000000000p-1, 0x1.419b423d5e800p-2, 0x1.8e436ec90e09dp-47}, - {0x1.7400000000000p-1, 0x1.4718dc271c800p-2, -0x1.f27ce0967d675p-45}, - {0x1.7200000000000p-1, 0x1.4c9e09e173000p-2, -0x1.e20891b0ad8a4p-45}, - {0x1.7000000000000p-1, 0x1.522ae0738a000p-2, 0x1.ebe708164c759p-45}, - {0x1.6e00000000000p-1, 0x1.57bf753c8d000p-2, 0x1.fadedee5d40efp-46}, - {0x1.6c00000000000p-1, 0x1.5d5bddf596000p-2, -0x1.a0b2a08a465dcp-47}, - }; - static const double A[] = { - -0x1p-1, - 0x1.555555555556p-2 * -2, - -0x1.0000000000006p-2 * -2, - 0x1.999999959554ep-3 * 4, - -0x1.555555529a47ap-3 * 4, - 0x1.2495b9b4845e9p-3 * -8, - -0x1.0002b8b263fc3p-3 * -8 - }; - static const double ln2hi = 0x1.62e42fefa3800p-1, - ln2lo = 0x1.ef35793c76730p-45; - - double z, r, y, invc, logc, logctail, kd, hi, t1, t2, lo, lo1, lo2, p; - double zhi, zlo, rhi, rlo, ar, ar2, ar3, lo3, lo4, arhi, arhi2; - UINT64 iz, tmp; - int k, i; - - /* x = 2^k z; where z is in range [OFF,2*OFF) and exact. - The range is split into N subintervals. - The ith subinterval contains z and c is near its center. */ - tmp = ix - 0x3fe6955500000000ULL; - i = (tmp >> (52 - 7)) % (1 << 7); - k = (INT64)tmp >> 52; /* arithmetic shift */ - iz = ix - (tmp & 0xfffULL << 52); - z = *(double*)&iz; - kd = k; - - /* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */ - invc = T[i].invc; - logc = T[i].logc; - logctail = T[i].logctail; - - /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and - |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */ - /* Split z such that rhi, rlo and rhi*rhi are exact and |rlo| <= |r|. */ - iz = (iz + (1ULL << 31)) & (-1ULL << 32); - zhi = *(double*)&iz; - zlo = z - zhi; - rhi = zhi * invc - 1.0; - rlo = zlo * invc; - r = rhi + rlo; - - /* k*Ln2 + log(c) + r. */ - t1 = kd * ln2hi + logc; - t2 = t1 + r; - lo1 = kd * ln2lo + logctail; - lo2 = t1 - t2 + r; - - /* Evaluation is optimized assuming superscalar pipelined execution. */ - ar = A[0] * r; /* A[0] = -0.5. */ - ar2 = r * ar; - ar3 = r * ar2; - /* k*Ln2 + log(c) + r + A[0]*r*r. */ - arhi = A[0] * rhi; - arhi2 = rhi * arhi; - hi = t2 + arhi2; - lo3 = rlo * (ar + arhi); - lo4 = t2 - hi + arhi2; - /* p = log1p(r) - r - A[0]*r*r. */ - p = (ar3 * (A[1] + r * A[2] + ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r * A[6])))); - lo = lo1 + lo2 + lo3 + lo4 + p; - y = hi + lo; - *tail = hi - y + lo; - return y; -} - -/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|. - The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1. */ -static double pow_exp(double argx, double argy, double x, double xtail, UINT32 sign_bias) -{ - static const double C[] = { - 0x1.ffffffffffdbdp-2, - 0x1.555555555543cp-3, - 0x1.55555cf172b91p-5, - 0x1.1111167a4d017p-7 - }; - static const double invln2N = 0x1.71547652b82fep0 * (1 << 7), - negln2hiN = -0x1.62e42fefa0000p-8, - negln2loN = -0x1.cf79abc9e3b3ap-47; - - UINT32 abstop; - UINT64 ki, idx, top, sbits; - double kd, z, r, r2, scale, tail, tmp; - - abstop = (*(UINT64*)&x >> 52) & 0x7ff; - if (abstop - 0x3c9 >= 0x408 - 0x3c9) { - if (abstop - 0x3c9 >= 0x80000000) { - /* Avoid spurious underflow for tiny x. */ - /* Note: 0 is common input. */ - double one = 1.0 + x; - return sign_bias ? -one : one; - } - if (abstop >= 0x409) { - /* Note: inf and nan are already handled. */ - if (*(UINT64*)&x >> 63) - return math_error(_UNDERFLOW, "pow", argx, argy, (sign_bias ? -DBL_MIN : DBL_MIN) * DBL_MIN); - return math_error(_OVERFLOW, "pow", argx, argy, (sign_bias ? -DBL_MAX : DBL_MAX) * DBL_MAX); - } - /* Large x is special cased below. */ - abstop = 0; - } - - /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */ - /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */ - z = invln2N * x; - kd = round(z); - ki = (INT64)kd; - r = x + kd * negln2hiN + kd * negln2loN; - /* The code assumes 2^-200 < |xtail| < 2^-8/N. */ - r += xtail; - /* 2^(k/N) ~= scale * (1 + tail). */ - idx = 2 * (ki % (1 << 7)); - top = (ki + sign_bias) << (52 - 7); - tail = *(double*)&exp_T[idx]; - /* This is only a valid scale when -1023*N < k < 1024*N. */ - sbits = exp_T[idx + 1] + top; - /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */ - /* Evaluation is optimized assuming superscalar pipelined execution. */ - r2 = r * r; - /* Without fma the worst case error is 0.25/N ulp larger. */ - /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp. */ - tmp = tail + r + r2 * (C[0] + r * C[1]) + r2 * r2 * (C[2] + r * C[3]); - if (abstop == 0) { - /* Handle cases that may overflow or underflow when computing the result that - is scale*(1+TMP) without intermediate rounding. The bit representation of - scale is in SBITS, however it has a computed exponent that may have - overflown into the sign bit so that needs to be adjusted before using it as - a double. (int32_t)KI is the k used in the argument reduction and exponent - adjustment of scale, positive k here means the result may overflow and - negative k means the result may underflow. */ - double scale, y; - - if ((ki & 0x80000000) == 0) { - /* k > 0, the exponent of scale might have overflowed by <= 460. */ - sbits -= 1009ull << 52; - scale = *(double*)&sbits; - y = 0x1p1009 * (scale + scale * tmp); - if (isinf(y)) - return math_error(_OVERFLOW, "pow", argx, argy, y); - return y; - } - /* k < 0, need special care in the subnormal range. */ - sbits += 1022ull << 52; - /* Note: sbits is signed scale. */ - scale = *(double*)&sbits; - y = scale + scale * tmp; - if (fabs(y) < 1.0) { - /* Round y to the right precision before scaling it into the subnormal - range to avoid double rounding that can cause 0.5+E/2 ulp error where - E is the worst-case ulp error outside the subnormal range. So this - is only useful if the goal is better than 1 ulp worst-case error. */ - double hi, lo, one = 1.0; - if (y < 0.0) - one = -1.0; - lo = scale - y + scale * tmp; - hi = one + y; - lo = one - hi + y + lo; - y = hi + lo - one; - /* Fix the sign of 0. */ - if (y == 0.0) { - sbits &= 0x8000000000000000ULL; - y = *(double*)&sbits; - } - /* The underflow exception needs to be signaled explicitly. */ - fp_barrier(fp_barrier(0x1p-1022) * 0x1p-1022); - y = 0x1p-1022 * y; - return math_error(_UNDERFLOW, "pow", argx, argy, y); - } - y = 0x1p-1022 * y; - return y; - } - scale = *(double*)&sbits; - /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there - is no spurious underflow here even without fma. */ - return scale + scale * tmp; -} - -/* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is - the bit representation of a non-zero finite floating-point value. */ -static inline int pow_checkint(UINT64 iy) -{ - int e = iy >> 52 & 0x7ff; - if (e < 0x3ff) - return 0; - if (e > 0x3ff + 52) - return 2; - if (iy & ((1ULL << (0x3ff + 52 - e)) - 1)) - return 0; - if (iy & (1ULL << (0x3ff + 52 - e))) - return 1; - return 2; -} - -/********************************************************************* - * pow (MSVCRT.@) - * - * Copied from musl: src/math/pow.c - */ -double CDECL pow( double x, double y ) -{ - UINT32 sign_bias = 0; - UINT64 ix, iy; - UINT32 topx, topy; - double lo, hi, ehi, elo, yhi, ylo, lhi, llo; - - ix = *(UINT64*)&x; - iy = *(UINT64*)&y; - topx = ix >> 52; - topy = iy >> 52; - if (topx - 0x001 >= 0x7ff - 0x001 || - (topy & 0x7ff) - 0x3be >= 0x43e - 0x3be) { - /* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0 - and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */ - /* Special cases: (x < 0x1p-126 or inf or nan) or - (|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */ - if (2 * iy - 1 >= 2 * 0x7ff0000000000000ULL - 1) { - if (2 * iy == 0) - return 1.0; - if (ix == 0x3ff0000000000000ULL) - return 1.0; - if (2 * ix > 2 * 0x7ff0000000000000ULL || - 2 * iy > 2 * 0x7ff0000000000000ULL) - return x + y; - if (2 * ix == 2 * 0x3ff0000000000000ULL) - return 1.0; - if ((2 * ix < 2 * 0x3ff0000000000000ULL) == !(iy >> 63)) - return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ - return y * y; - } - if (2 * ix - 1 >= 2 * 0x7ff0000000000000ULL - 1) { - double x2 = x * x; - if (ix >> 63 && pow_checkint(iy) == 1) - x2 = -x2; - if (iy & 0x8000000000000000ULL && x2 == 0.0) - return math_error(_SING, "pow", x, y, 1 / x2); - /* Without the barrier some versions of clang hoist the 1/x2 and - thus division by zero exception can be signaled spuriously. */ - return iy >> 63 ? fp_barrier(1 / x2) : x2; - } - /* Here x and y are non-zero finite. */ - if (ix >> 63) { - /* Finite x < 0. */ - int yint = pow_checkint(iy); - if (yint == 0) - return math_error(_DOMAIN, "pow", x, y, 0 / (x - x)); - if (yint == 1) - sign_bias = 0x800 << 7; - ix &= 0x7fffffffffffffff; - topx &= 0x7ff; - } - if ((topy & 0x7ff) - 0x3be >= 0x43e - 0x3be) { - /* Note: sign_bias == 0 here because y is not odd. */ - if (ix == 0x3ff0000000000000ULL) - return 1.0; - if ((topy & 0x7ff) < 0x3be) { - /* |y| < 2^-65, x^y ~= 1 + y*log(x). */ - return ix > 0x3ff0000000000000ULL ? 1.0 + y : 1.0 - y; - } - if ((ix > 0x3ff0000000000000ULL) == (topy < 0x800)) - return math_error(_OVERFLOW, "pow", x, y, fp_barrier(DBL_MAX) * DBL_MAX); - return math_error(_UNDERFLOW, "pow", x, y, fp_barrier(DBL_MIN) * DBL_MIN); - } - if (topx == 0) { - /* Normalize subnormal x so exponent becomes negative. */ - x *= 0x1p52; - ix = *(UINT64*)&x; - ix &= 0x7fffffffffffffff; - ix -= 52ULL << 52; - } - } - - hi = pow_log(ix, &lo); - iy &= -1ULL << 27; - yhi = *(double*)&iy; - ylo = y - yhi; - *(UINT64*)&lhi = *(UINT64*)&hi & -1ULL << 27; - llo = fp_barrier(hi - lhi + lo); - ehi = yhi * lhi; - elo = ylo * lhi + y * llo; /* |elo| < |ehi| * 2^-25. */ - return pow_exp(x, y, ehi, elo, sign_bias); -} - static BOOL sqrt_validate( double *x, BOOL update_sw ) { short c = _dclass(*x); diff --git a/libs/musl/src/math/pow.c b/libs/musl/src/math/pow.c index cce83687eb9..3b7815c1317 100644 --- a/libs/musl/src/math/pow.c +++ b/libs/musl/src/math/pow.c @@ -121,7 +121,7 @@ static inline double_t log_inline(uint64_t ix, double_t *tail) a double. (int32_t)KI is the k used in the argument reduction and exponent adjustment of scale, positive k here means the result may overflow and negative k means the result may underflow. */ -static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki) +static inline double specialcase(double argx, double argy, double_t tmp, uint64_t sbits, uint64_t ki) { double_t scale, y;
@@ -130,6 +130,8 @@ static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki) sbits -= 1009ull << 52; scale = asdouble(sbits); y = 0x1p1009 * (scale + scale * tmp); + if (isinf(y)) + return math_error(_OVERFLOW, "pow", argx, argy, y); return eval_as_double(y); } /* k < 0, need special care in the subnormal range. */ @@ -154,6 +156,8 @@ static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki) y = asdouble(sbits & 0x8000000000000000); /* The underflow exception needs to be signaled explicitly. */ fp_force_eval(fp_barrier(0x1p-1022) * 0x1p-1022); + y = 0x1p-1022 * y; + return math_error(_UNDERFLOW, "pow", argx, argy, y); } y = 0x1p-1022 * y; return eval_as_double(y); @@ -163,7 +167,7 @@ static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki)
/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|. The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1. */ -static inline double exp_inline(double_t x, double_t xtail, uint32_t sign_bias) +static inline double exp_inline(double argx, double argy, double_t x, double_t xtail, uint32_t sign_bias) { uint32_t abstop; uint64_t ki, idx, top, sbits; @@ -182,9 +186,9 @@ static inline double exp_inline(double_t x, double_t xtail, uint32_t sign_bias) if (abstop >= top12(1024.0)) { /* Note: inf and nan are already handled. */ if (asuint64(x) >> 63) - return __math_uflow(sign_bias); + return math_error(_UNDERFLOW, "pow", argx, argy, (sign_bias ? -DBL_MIN : DBL_MIN) * DBL_MIN); else - return __math_oflow(sign_bias); + return math_error(_OVERFLOW, "pow", argx, argy, (sign_bias ? -DBL_MAX : DBL_MAX) * DBL_MAX); } /* Large x is special cased below. */ abstop = 0; @@ -193,20 +197,8 @@ static inline double exp_inline(double_t x, double_t xtail, uint32_t sign_bias) /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */ /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */ z = InvLn2N * x; -#if TOINT_INTRINSICS - kd = roundtoint(z); - ki = converttoint(z); -#elif EXP_USE_TOINT_NARROW - /* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes. */ - kd = eval_as_double(z + Shift); - ki = asuint64(kd) >> 16; - kd = (double_t)(int32_t)ki; -#else - /* z - kd is in [-1, 1] in non-nearest rounding modes. */ - kd = eval_as_double(z + Shift); - ki = asuint64(kd); - kd -= Shift; -#endif + kd = round(z); + ki = (int64_t)kd; r = x + kd * NegLn2hiN + kd * NegLn2loN; /* The code assumes 2^-200 < |xtail| < 2^-8/N. */ r += xtail; @@ -223,7 +215,7 @@ static inline double exp_inline(double_t x, double_t xtail, uint32_t sign_bias) /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp. */ tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5); if (predict_false(abstop == 0)) - return specialcase(tmp, sbits, ki); + return specialcase(argx, argy, tmp, sbits, ki); scale = asdouble(sbits); /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there is no spurious underflow here even without fma. */ @@ -286,6 +278,8 @@ double __cdecl pow(double x, double y) double_t x2 = x * x; if (ix >> 63 && checkint(iy) == 1) x2 = -x2; + if (iy & 0x8000000000000000ULL && x2 == 0.0) + return math_error(_SING, "pow", x, y, 1 / x2); /* Without the barrier some versions of clang hoist the 1/x2 and thus division by zero exception can be signaled spuriously. */ return iy >> 63 ? fp_barrier(1 / x2) : x2; @@ -295,7 +289,7 @@ double __cdecl pow(double x, double y) /* Finite x < 0. */ int yint = checkint(iy); if (yint == 0) - return __math_invalid(x); + return math_error(_DOMAIN, "pow", x, y, 0 / (x - x)); if (yint == 1) sign_bias = SIGN_BIAS; ix &= 0x7fffffffffffffff; @@ -313,9 +307,9 @@ double __cdecl pow(double x, double y) else return 1.0; } - return (ix > asuint64(1.0)) == (topy < 0x800) ? - __math_oflow(0) : - __math_uflow(0); + if ((ix > asuint64(1.0)) == (topy < 0x800)) + return math_error(_OVERFLOW, "pow", x, y, fp_barrier(DBL_MAX) * DBL_MAX); + return math_error(_UNDERFLOW, "pow", x, y, fp_barrier(DBL_MIN) * DBL_MIN); } if (topx == 0) { /* Normalize subnormal x so exponent becomes negative. */ @@ -339,5 +333,5 @@ double __cdecl pow(double x, double y) ehi = yhi * lhi; elo = ylo * lhi + y * llo; /* |elo| < |ehi| * 2^-25. */ #endif - return exp_inline(ehi, elo, sign_bias); + return exp_inline(x, y, ehi, elo, sign_bias); } diff --git a/libs/musl/src/math/powf.c b/libs/musl/src/math/powf.c index fba6270713b..2e8576ba49f 100644 --- a/libs/musl/src/math/powf.c +++ b/libs/musl/src/math/powf.c @@ -73,19 +73,10 @@ static inline float exp2_inline(double_t xd, uint32_t sign_bias) uint64_t ki, ski, t; double_t kd, z, r, r2, y, s;
-#if TOINT_INTRINSICS #define C __exp2f_data.poly_scaled /* N*x = k + r with r in [-1/2, 1/2] */ - kd = roundtoint(xd); /* k */ - ki = converttoint(xd); -#else -#define C __exp2f_data.poly -#define SHIFT __exp2f_data.shift_scaled - /* x = k/N + r with r in [-1/(2N), 1/(2N)] */ - kd = eval_as_double(xd + SHIFT); - ki = asuint64(kd); - kd -= SHIFT; /* k/N */ -#endif + kd = round(xd); /* k */ + ki = (int64_t)kd; r = xd - kd;
/* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ @@ -150,6 +141,8 @@ float __cdecl powf(float x, float y) float_t x2 = x * x; if (ix & 0x80000000 && checkint(iy) == 1) x2 = -x2; + if (iy & 0x80000000 && x2 == 0.0) + return math_error(_SING, "powf", x, y, 1 / x2); /* Without the barrier some versions of clang hoist the 1/x2 and thus division by zero exception can be signaled spuriously. */ return iy & 0x80000000 ? fp_barrierf(1 / x2) : x2; @@ -159,7 +152,7 @@ float __cdecl powf(float x, float y) /* Finite x < 0. */ int yint = checkint(iy); if (yint == 0) - return __math_invalidf(x); + return math_error(_DOMAIN, "powf", x, y, 0 / (x - x)); if (yint == 1) sign_bias = SIGN_BIAS; ix &= 0x7fffffff; @@ -177,9 +170,9 @@ float __cdecl powf(float x, float y) asuint64(126.0 * POWF_SCALE) >> 47)) { /* |y*log(x)| >= 126. */ if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE) - return __math_oflowf(sign_bias); + return math_error(_OVERFLOW, "powf", x, y, (sign_bias ? -1.0 : 1.0) * 0x1p1023); if (ylogx <= -150.0 * POWF_SCALE) - return __math_uflowf(sign_bias); + return math_error(_UNDERFLOW, "powf", x, y, (sign_bias ? -1.0 : 1.0) / 0x1p1023); } return exp2_inline(ylogx, sign_bias); }