diff options
Diffstat (limited to 'Modules/mathmodule.c')
| -rw-r--r-- | Modules/mathmodule.c | 25 |
1 files changed, 22 insertions, 3 deletions
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 11d9b7418a2..71d9c1387f5 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -2008,13 +2008,11 @@ math.factorial / Find n!. - -Raise a ValueError if x is negative or non-integral. [clinic start generated code]*/ static PyObject * math_factorial(PyObject *module, PyObject *arg) -/*[clinic end generated code: output=6686f26fae00e9ca input=713fb771677e8c31]*/ +/*[clinic end generated code: output=6686f26fae00e9ca input=366cc321df3d4773]*/ { long x, two_valuation; int overflow; @@ -2163,6 +2161,27 @@ math_ldexp_impl(PyObject *module, double x, PyObject *i) } else { errno = 0; r = ldexp(x, (int)exp); +#ifdef _MSC_VER + if (DBL_MIN > r && r > -DBL_MIN) { + /* Denormal (or zero) results can be incorrectly rounded here (rather, + truncated). Fixed in newer versions of the C runtime, included + with Windows 11. */ + int original_exp; + frexp(x, &original_exp); + if (original_exp > DBL_MIN_EXP) { + /* Shift down to the smallest normal binade. No bits lost. */ + int shift = DBL_MIN_EXP - original_exp; + x = ldexp(x, shift); + exp -= shift; + } + /* Multiplying by 2**exp finishes the job, and the HW will round as + appropriate. Note: if exp < -DBL_MANT_DIG, all of x is shifted + to be < 0.5ULP of smallest denorm, so should be thrown away. If + exp is so very negative that ldexp underflows to 0, that's fine; + no need to check in advance. */ + r = x*ldexp(1.0, (int)exp); + } +#endif if (isinf(r)) errno = ERANGE; } |