Use fused multiply add to improve logarithms of bigints

[tim-one]

Feature or enhancement

Proposal:

This line from mathmodule.c's loghelper():

            result = func(x) + func(2.0) * e;

is begging to be written to use a HW fma instruction. There isn't massive cancellation here, where fma can "do magic", just the ordinary benefit of suffering one rounding error instead of 2. A quick Python prototype suggests it would do real good.

Here are summaries of two runs.

• Testing log10() on "random" bigints of up to 100 million decimal digits.
• Both runs use the same inputs.
• The result of mpmath.log10(), with prec=100, is taken to be correct.
• A line starting with, say, 0.4, records the number of cases where the computed result was <= 0.4 ULP too large, but greater than 0.3 ULP too large.

Before (CPython 3.14):
-0.8546426777998732 ULP to 0.9945997134441882 ULP
-0.8    9  0.2% **
-0.7   52  1.0% *******
-0.6  118  2.4% **************
-0.5  233  4.7% ****************************
-0.4  261  5.2% *******************************
-0.3  299  6.0% ***********************************
-0.2  318  6.4% **************************************
-0.1  348  7.0% *****************************************
 0.0  480  9.6% *********************************************************
 0.1  514 10.3% ************************************************************
 0.2  506 10.1% ************************************************************
 0.3  419  8.4% *************************************************
 0.4  396  7.9% ***********************************************
 0.5  326  6.5% ***************************************
 0.6  267  5.3% ********************************
 0.7  171  3.4% ********************
 0.8  132  2.6% ****************
 0.9  118  2.4% **************
 1.0   33  0.7% ****

After:
-0.454879880773845 ULP to 0.5800814567479762 ULP
-0.4  234  4.7% **************************
-0.3  470  9.4% ***************************************************
-0.2  450  9.0% *************************************************
-0.1  466  9.3% ***************************************************
 0.0  513 10.3% ********************************************************
 0.1  514 10.3% ********************************************************
 0.2  515 10.3% ********************************************************
 0.3  471  9.4% ***************************************************
 0.4  514 10.3% ********************************************************
 0.5  559 11.2% ************************************************************
 0.6  294  5.9% ********************************

So max errors (in both directions) were cut, and the error distribution got significantly tighter.

Note that, either way, our results tend to be "too large". I believe that's because, on my Windows box, while log10(2) returns the best possible (representable) result, it's nevertheless larger than the infinitely precise value, and multiplying by e amplifies that.

This is a very easy one if someone wants to beat me to it, although, ya, testing can consume as much of your future life as you allow it to 😉.

Has this already been discussed elsewhere?

This is a minor feature, which does not need previous discussion elsewhere

Links to previous discussion of this feature:

No response

Linked PRs

• gh-140443: Use fma in loghelper to improve accuracy of log for very large integers #140469

[Abhi210]

Hello! If @tim-one does not mind, I would love to take over this issue. Thank You!

[tim-one]

Be my guest! I think it's a great way to ease into CPython development (small, well-defined, conceptually clear). Have fun 😄. Note: in the C code, you'll have to cast e to double to prevent compiler warnings.

[serhiy-storchaka]

Does it work the same for much smaller ranges?

[tim-one]

Sorry, I don't understand your question. But maybe this will answer it: the code in question is invoked if and only if converting an int to a double overflows (in which case the overflow is suppressed, and we take an entirely path specific to "very large" ints, like over (roughly) 2**1024).

[serhiy-storchaka]

Thank you for explanation. I forgot about this.

[Abhi210]

Hello, so I changed the line in loghelper to use fma according to:

result = fma(func(2.0),(double)e,func(x));

For testing, I reproduced the ULP error distribution using a Python script modeled on the one you described earlier. The main bottleneck is the mpmath.log10 call at mpmath.mp.prec = 100, which makes the test run quite slowly for large random integers.

Here’s the test I used:

def test_ulp_distribution(num_tests=1000, max_digits=100000):

    random.seed(12345)

    import mpmath
    mpmath.mp.prec = 100  # 100 decimal digits precision

    ulp_errors = []

    for i in range(num_tests):
        # Generate random bigint
        #num_digits = random.randint(100, max_digits)
        #num = random.randint(10**(num_digits-1), 10**num_digits - 1)

        # Force integers too large for float conversion
        num_bits = random.randint(1024, 10_000_000)  # 1024 bits ≈ 308 digits
        mantissa = random.randint(1 << 52, (1 << 53) - 1)  # 53-bit mantissa
        num = mantissa << (num_bits - 53)

        # Get Python's result
        py_result = math.log10(num)

        # Get high-precision reference
        mp_result = float(mpmath.log10(num))

        # Calculate ULP error
        if math.isfinite(mp_result):
            ulp_error = (py_result - mp_result) / math.ulp(mp_result)
            ulp_errors.append(ulp_error)

    # Create histogram
    print(f"\nULP Error Range: {min(ulp_errors):.2f} to {max(ulp_errors):.2f}")

    # Create bins
    bins = {}
    for ulp in ulp_errors:
        bin_key = round(ulp * 10) / 10 # Round to nearest 0.1
        bins[bin_key] = bins.get(bin_key, 0) + 1

    # Print histogram
    print("\nULP Error Distribution:")
    for bin_key in sorted(bins.keys()):
        count = bins[bin_key]
        pct = 100.0 * count / len(ulp_errors)
        bar = '*' * int(pct * 2)  # Scale for display
        print(f"{bin_key:6.2f} {count:6d} {pct:6.2f}% {bar}")

It’s currently running (about halfway through) and seems stable so far.

Could you please confirm if this approach looks correct to you?
Also, are there any additional cases or precision settings you’d recommend testing (for example, smaller mp.prec or specific edge values)?

Thanks!

[StanFromIreland]

result = fma(func(2.0),(double)e,func(x));

Looks correct to me, although I'm not Tim:-) Please do include spaces:

- result = fma(func(2.0),(double)e,func(x));
+ result = fma(func(2.0), (double)e, func(x));