ntdll: Copy pow() implementation from msvcrt.
Signed-off-by: Alexandre Julliard <julliard@winehq.org>
This commit is contained in:
Alexandre Julliard
parent
0c0915e985
commit
140cd3553c
|
|
@ -315,6 +315,27 @@ recompute:
|
|||
return n & 7;
|
||||
}
|
||||
|
||||
/* Based on musl implementation: src/math/round.c */
|
||||
static double __round(double x)
|
||||
{
|
||||
ULONGLONG llx = *(ULONGLONG*)&x, tmp;
|
||||
int e = (llx >> 52 & 0x7ff) - 0x3ff;
|
||||
|
||||
if (e >= 52)
|
||||
return x;
|
||||
if (e < -1)
|
||||
return 0 * x;
|
||||
else if (e == -1)
|
||||
return signbit(x) ? -1 : 1;
|
||||
|
||||
tmp = 0x000fffffffffffffULL >> e;
|
||||
if (!(llx & tmp))
|
||||
return x;
|
||||
llx += 0x0008000000000000ULL >> e;
|
||||
llx &= ~tmp;
|
||||
return *(double*)&llx;
|
||||
}
|
||||
|
||||
/* Copied from musl: src/math/__rem_pio2.c */
|
||||
static int __rem_pio2(double x, double *y)
|
||||
{
|
||||
|
|
@ -566,6 +587,466 @@ static double __tan(double x, double y, int odd)
|
|||
return a0 + a * (1.0 + a0 * w0 + a0 * v);
|
||||
}
|
||||
|
||||
/* Copied from musl: src/math/exp_data.c */
|
||||
static const UINT64 exp_T[] = {
|
||||
0x0ULL, 0x3ff0000000000000ULL,
|
||||
0x3c9b3b4f1a88bf6eULL, 0x3feff63da9fb3335ULL,
|
||||
0xbc7160139cd8dc5dULL, 0x3fefec9a3e778061ULL,
|
||||
0xbc905e7a108766d1ULL, 0x3fefe315e86e7f85ULL,
|
||||
0x3c8cd2523567f613ULL, 0x3fefd9b0d3158574ULL,
|
||||
0xbc8bce8023f98efaULL, 0x3fefd06b29ddf6deULL,
|
||||
0x3c60f74e61e6c861ULL, 0x3fefc74518759bc8ULL,
|
||||
0x3c90a3e45b33d399ULL, 0x3fefbe3ecac6f383ULL,
|
||||
0x3c979aa65d837b6dULL, 0x3fefb5586cf9890fULL,
|
||||
0x3c8eb51a92fdeffcULL, 0x3fefac922b7247f7ULL,
|
||||
0x3c3ebe3d702f9cd1ULL, 0x3fefa3ec32d3d1a2ULL,
|
||||
0xbc6a033489906e0bULL, 0x3fef9b66affed31bULL,
|
||||
0xbc9556522a2fbd0eULL, 0x3fef9301d0125b51ULL,
|
||||
0xbc5080ef8c4eea55ULL, 0x3fef8abdc06c31ccULL,
|
||||
0xbc91c923b9d5f416ULL, 0x3fef829aaea92de0ULL,
|
||||
0x3c80d3e3e95c55afULL, 0x3fef7a98c8a58e51ULL,
|
||||
0xbc801b15eaa59348ULL, 0x3fef72b83c7d517bULL,
|
||||
0xbc8f1ff055de323dULL, 0x3fef6af9388c8deaULL,
|
||||
0x3c8b898c3f1353bfULL, 0x3fef635beb6fcb75ULL,
|
||||
0xbc96d99c7611eb26ULL, 0x3fef5be084045cd4ULL,
|
||||
0x3c9aecf73e3a2f60ULL, 0x3fef54873168b9aaULL,
|
||||
0xbc8fe782cb86389dULL, 0x3fef4d5022fcd91dULL,
|
||||
0x3c8a6f4144a6c38dULL, 0x3fef463b88628cd6ULL,
|
||||
0x3c807a05b0e4047dULL, 0x3fef3f49917ddc96ULL,
|
||||
0x3c968efde3a8a894ULL, 0x3fef387a6e756238ULL,
|
||||
0x3c875e18f274487dULL, 0x3fef31ce4fb2a63fULL,
|
||||
0x3c80472b981fe7f2ULL, 0x3fef2b4565e27cddULL,
|
||||
0xbc96b87b3f71085eULL, 0x3fef24dfe1f56381ULL,
|
||||
0x3c82f7e16d09ab31ULL, 0x3fef1e9df51fdee1ULL,
|
||||
0xbc3d219b1a6fbffaULL, 0x3fef187fd0dad990ULL,
|
||||
0x3c8b3782720c0ab4ULL, 0x3fef1285a6e4030bULL,
|
||||
0x3c6e149289cecb8fULL, 0x3fef0cafa93e2f56ULL,
|
||||
0x3c834d754db0abb6ULL, 0x3fef06fe0a31b715ULL,
|
||||
0x3c864201e2ac744cULL, 0x3fef0170fc4cd831ULL,
|
||||
0x3c8fdd395dd3f84aULL, 0x3feefc08b26416ffULL,
|
||||
0xbc86a3803b8e5b04ULL, 0x3feef6c55f929ff1ULL,
|
||||
0xbc924aedcc4b5068ULL, 0x3feef1a7373aa9cbULL,
|
||||
0xbc9907f81b512d8eULL, 0x3feeecae6d05d866ULL,
|
||||
0xbc71d1e83e9436d2ULL, 0x3feee7db34e59ff7ULL,
|
||||
0xbc991919b3ce1b15ULL, 0x3feee32dc313a8e5ULL,
|
||||
0x3c859f48a72a4c6dULL, 0x3feedea64c123422ULL,
|
||||
0xbc9312607a28698aULL, 0x3feeda4504ac801cULL,
|
||||
0xbc58a78f4817895bULL, 0x3feed60a21f72e2aULL,
|
||||
0xbc7c2c9b67499a1bULL, 0x3feed1f5d950a897ULL,
|
||||
0x3c4363ed60c2ac11ULL, 0x3feece086061892dULL,
|
||||
0x3c9666093b0664efULL, 0x3feeca41ed1d0057ULL,
|
||||
0x3c6ecce1daa10379ULL, 0x3feec6a2b5c13cd0ULL,
|
||||
0x3c93ff8e3f0f1230ULL, 0x3feec32af0d7d3deULL,
|
||||
0x3c7690cebb7aafb0ULL, 0x3feebfdad5362a27ULL,
|
||||
0x3c931dbdeb54e077ULL, 0x3feebcb299fddd0dULL,
|
||||
0xbc8f94340071a38eULL, 0x3feeb9b2769d2ca7ULL,
|
||||
0xbc87deccdc93a349ULL, 0x3feeb6daa2cf6642ULL,
|
||||
0xbc78dec6bd0f385fULL, 0x3feeb42b569d4f82ULL,
|
||||
0xbc861246ec7b5cf6ULL, 0x3feeb1a4ca5d920fULL,
|
||||
0x3c93350518fdd78eULL, 0x3feeaf4736b527daULL,
|
||||
0x3c7b98b72f8a9b05ULL, 0x3feead12d497c7fdULL,
|
||||
0x3c9063e1e21c5409ULL, 0x3feeab07dd485429ULL,
|
||||
0x3c34c7855019c6eaULL, 0x3feea9268a5946b7ULL,
|
||||
0x3c9432e62b64c035ULL, 0x3feea76f15ad2148ULL,
|
||||
0xbc8ce44a6199769fULL, 0x3feea5e1b976dc09ULL,
|
||||
0xbc8c33c53bef4da8ULL, 0x3feea47eb03a5585ULL,
|
||||
0xbc845378892be9aeULL, 0x3feea34634ccc320ULL,
|
||||
0xbc93cedd78565858ULL, 0x3feea23882552225ULL,
|
||||
0x3c5710aa807e1964ULL, 0x3feea155d44ca973ULL,
|
||||
0xbc93b3efbf5e2228ULL, 0x3feea09e667f3bcdULL,
|
||||
0xbc6a12ad8734b982ULL, 0x3feea012750bdabfULL,
|
||||
0xbc6367efb86da9eeULL, 0x3fee9fb23c651a2fULL,
|
||||
0xbc80dc3d54e08851ULL, 0x3fee9f7df9519484ULL,
|
||||
0xbc781f647e5a3ecfULL, 0x3fee9f75e8ec5f74ULL,
|
||||
0xbc86ee4ac08b7db0ULL, 0x3fee9f9a48a58174ULL,
|
||||
0xbc8619321e55e68aULL, 0x3fee9feb564267c9ULL,
|
||||
0x3c909ccb5e09d4d3ULL, 0x3feea0694fde5d3fULL,
|
||||
0xbc7b32dcb94da51dULL, 0x3feea11473eb0187ULL,
|
||||
0x3c94ecfd5467c06bULL, 0x3feea1ed0130c132ULL,
|
||||
0x3c65ebe1abd66c55ULL, 0x3feea2f336cf4e62ULL,
|
||||
0xbc88a1c52fb3cf42ULL, 0x3feea427543e1a12ULL,
|
||||
0xbc9369b6f13b3734ULL, 0x3feea589994cce13ULL,
|
||||
0xbc805e843a19ff1eULL, 0x3feea71a4623c7adULL,
|
||||
0xbc94d450d872576eULL, 0x3feea8d99b4492edULL,
|
||||
0x3c90ad675b0e8a00ULL, 0x3feeaac7d98a6699ULL,
|
||||
0x3c8db72fc1f0eab4ULL, 0x3feeace5422aa0dbULL,
|
||||
0xbc65b6609cc5e7ffULL, 0x3feeaf3216b5448cULL,
|
||||
0x3c7bf68359f35f44ULL, 0x3feeb1ae99157736ULL,
|
||||
0xbc93091fa71e3d83ULL, 0x3feeb45b0b91ffc6ULL,
|
||||
0xbc5da9b88b6c1e29ULL, 0x3feeb737b0cdc5e5ULL,
|
||||
0xbc6c23f97c90b959ULL, 0x3feeba44cbc8520fULL,
|
||||
0xbc92434322f4f9aaULL, 0x3feebd829fde4e50ULL,
|
||||
0xbc85ca6cd7668e4bULL, 0x3feec0f170ca07baULL,
|
||||
0x3c71affc2b91ce27ULL, 0x3feec49182a3f090ULL,
|
||||
0x3c6dd235e10a73bbULL, 0x3feec86319e32323ULL,
|
||||
0xbc87c50422622263ULL, 0x3feecc667b5de565ULL,
|
||||
0x3c8b1c86e3e231d5ULL, 0x3feed09bec4a2d33ULL,
|
||||
0xbc91bbd1d3bcbb15ULL, 0x3feed503b23e255dULL,
|
||||
0x3c90cc319cee31d2ULL, 0x3feed99e1330b358ULL,
|
||||
0x3c8469846e735ab3ULL, 0x3feede6b5579fdbfULL,
|
||||
0xbc82dfcd978e9db4ULL, 0x3feee36bbfd3f37aULL,
|
||||
0x3c8c1a7792cb3387ULL, 0x3feee89f995ad3adULL,
|
||||
0xbc907b8f4ad1d9faULL, 0x3feeee07298db666ULL,
|
||||
0xbc55c3d956dcaebaULL, 0x3feef3a2b84f15fbULL,
|
||||
0xbc90a40e3da6f640ULL, 0x3feef9728de5593aULL,
|
||||
0xbc68d6f438ad9334ULL, 0x3feeff76f2fb5e47ULL,
|
||||
0xbc91eee26b588a35ULL, 0x3fef05b030a1064aULL,
|
||||
0x3c74ffd70a5fddcdULL, 0x3fef0c1e904bc1d2ULL,
|
||||
0xbc91bdfbfa9298acULL, 0x3fef12c25bd71e09ULL,
|
||||
0x3c736eae30af0cb3ULL, 0x3fef199bdd85529cULL,
|
||||
0x3c8ee3325c9ffd94ULL, 0x3fef20ab5fffd07aULL,
|
||||
0x3c84e08fd10959acULL, 0x3fef27f12e57d14bULL,
|
||||
0x3c63cdaf384e1a67ULL, 0x3fef2f6d9406e7b5ULL,
|
||||
0x3c676b2c6c921968ULL, 0x3fef3720dcef9069ULL,
|
||||
0xbc808a1883ccb5d2ULL, 0x3fef3f0b555dc3faULL,
|
||||
0xbc8fad5d3ffffa6fULL, 0x3fef472d4a07897cULL,
|
||||
0xbc900dae3875a949ULL, 0x3fef4f87080d89f2ULL,
|
||||
0x3c74a385a63d07a7ULL, 0x3fef5818dcfba487ULL,
|
||||
0xbc82919e2040220fULL, 0x3fef60e316c98398ULL,
|
||||
0x3c8e5a50d5c192acULL, 0x3fef69e603db3285ULL,
|
||||
0x3c843a59ac016b4bULL, 0x3fef7321f301b460ULL,
|
||||
0xbc82d52107b43e1fULL, 0x3fef7c97337b9b5fULL,
|
||||
0xbc892ab93b470dc9ULL, 0x3fef864614f5a129ULL,
|
||||
0x3c74b604603a88d3ULL, 0x3fef902ee78b3ff6ULL,
|
||||
0x3c83c5ec519d7271ULL, 0x3fef9a51fbc74c83ULL,
|
||||
0xbc8ff7128fd391f0ULL, 0x3fefa4afa2a490daULL,
|
||||
0xbc8dae98e223747dULL, 0x3fefaf482d8e67f1ULL,
|
||||
0x3c8ec3bc41aa2008ULL, 0x3fefba1bee615a27ULL,
|
||||
0x3c842b94c3a9eb32ULL, 0x3fefc52b376bba97ULL,
|
||||
0x3c8a64a931d185eeULL, 0x3fefd0765b6e4540ULL,
|
||||
0xbc8e37bae43be3edULL, 0x3fefdbfdad9cbe14ULL,
|
||||
0x3c77893b4d91cd9dULL, 0x3fefe7c1819e90d8ULL,
|
||||
0x3c5305c14160cc89ULL, 0x3feff3c22b8f71f1ULL
|
||||
};
|
||||
|
||||
/* 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 (sign_bias ? -DBL_MIN : DBL_MIN) * DBL_MIN;
|
||||
return (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);
|
||||
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 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;
|
||||
}
|
||||
|
||||
|
||||
/*********************************************************************
|
||||
* abs (NTDLL.@)
|
||||
|
|
@ -1146,10 +1627,91 @@ double CDECL log( double x )
|
|||
|
||||
/*********************************************************************
|
||||
* pow (NTDLL.@)
|
||||
*
|
||||
* Copied from musl: src/math/pow.c
|
||||
*/
|
||||
double CDECL pow( double x, double y )
|
||||
{
|
||||
return unix_funcs->pow( x, 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 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 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 fp_barrier(DBL_MAX) * DBL_MAX;
|
||||
return 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);
|
||||
}
|
||||
|
||||
/*********************************************************************
|
||||
|
|
|
|||
Loading…
Reference in New Issue