gh-140443: Use `fma` in `loghelper` to improve accuracy of log for very large integer...

* gh-140443:use fma in loghelper to improve accuracy of log for very large integers

Use fused multiply-add in log_helper() for huge ints.

Saving a rounding here is remarkably effective. Across some millions
of randomized test cases with ints up to a billion bits, on Windows
and using log10, the ULP error distribution was dramatically
flattened, and its range was nearly cut in half. In fact, the largest
error Tim saw was under 0.6 ULP.

---------

Co-authored-by: abhi210 <27881020+Abhi210@users.noreply.github.com>
Co-authored-by: Stan Ulbrych <89152624+StanFromIreland@users.noreply.github.com>

diff --git a/Misc/ACKS b/Misc/ACKS
index 6876380e0ba8d2097f2d3484511ecfb1930a286c..f5f15f2eb7ea24341432955e94674ea56f38aad1 100644 (file)
--- a/Misc/ACKS
+++ b/Misc/ACKS
@@ -1921,6 +1921,7 @@ Tim Tisdall
 Jason Tishler
 Christian Tismer
 Jim Tittsler
+Abhishek Tiwari
 Frank J. Tobin
 James Tocknell
 Bennett Todd

diff --git a/Misc/NEWS.d/next/Core_and_Builtins/2025-10-22-23-26-37.gh-issue-140443.wT5i1A.rst b/Misc/NEWS.d/next/Core_and_Builtins/2025-10-22-23-26-37.gh-issue-140443.wT5i1A.rst
new file mode 100644 (file)
index 0000000..a1fff8f
--- /dev/null
+++ b/Misc/NEWS.d/next/Core_and_Builtins/2025-10-22-23-26-37.gh-issue-140443.wT5i1A.rst
@@ -0,0 +1,5 @@
+The logarithm functions (such as :func:`math.log10` and :func:`math.log`) may now produce
+slightly different results for extremely large integers that cannot be
+converted to floats without overflow. These results are generally more
+accurate, with reduced worst-case error and a tighter overall error
+distribution.

diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index c631beb9ce5477fafc5bc242b19456f3658a7042..be88841716b00416b3072d3078eafb13394f2c97 100644 (file)
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2309,7 +2309,7 @@ loghelper(PyObject* arg, double (*func)(double))
             assert(e >= 0);
             assert(!PyErr_Occurred());
             /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
-            result = func(x) + func(2.0) * e;
+            result = fma(func(2.0), (double)e, func(x));
         }
         else
             /* Successfully converted x to a double. */