gh-132876: workaround broken ldexp() on Windows 10

Co-authored-by: Tim Peters <[email protected]>

Author: skirpichev

Misc/NEWS.d/next/Library/2025-04-29-11-48-46.gh-issue-132876.lyTQGZ.rst

@@ -0,0 +1 @@
+Workaround broken ``ldexp()`` implementation on Windows 10.

Modules/mathmodule.c

@@ -2163,6 +2163,26 @@ math_ldexp_impl(PyObject *module, double x, PyObject *i)
     } else {
         errno = 0;
         r = ldexp(x, (int)exp);
+#if _MSC_VER
+        if (r && r > -DBL_MIN && r < DBL_MIN) {
+            /* Denormal result 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. */
+                x = ldexp(x, DBL_MIN_EXP - original_exp);
+                exp += original_exp - DBL_MIN_EXP;
+            }
+            /* 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;
     }