math.ldexp gives incorrect results on Windows

[tkreuzer]

Bug report

Bug description:

>>> import math
>>> math.ldexp(6993274598585239, -1126)
5e-324
>>>

The correct result would be 1e-323. This is obviously a bug in the Windows ldexp implementation (it works fine on Linux). But it would be good, if it could be fixed on the CPython side with a workaround.
math.ldexp is used by mpmath to round from multiprecision to float, leading to incorrect rounding on Windows.

CPython versions tested on:

3.13

Operating systems tested on:

Windows

Linked PRs

• gh-132876: workaround broken ldexp() on Windows 10 #133135
• [3.14] gh-132876: workaround broken ldexp() on Windows 10 (GH-133135) #134684
• [3.13] gh-132876: workaround broken ldexp() on Windows 10 (GH-133135) #134685

[skirpichev]

For a context: mpmath/mpmath#942

[mmmp-oo]

Hi! 👋

This issue appears to be caused by platform-specific differences in floating-point handling, particularly on Windows. For consistency across platforms, a workaround is to use arbitrary-precision libraries such as mpmath:

from mpmath import mp
mp.dps = 50
x = mp.mpf(6993274598585239)
y = -1126
print(mp.ldexp(x, y))

Of course, since math.ldexp is a native function, a permanent solution would need to be implemented within CPython itself.

Hope this helps! 😊

[skirpichev]

appears to be caused by platform-specific differences

Just a bug in the platform libc. @tkreuzer, could you please describe your system (OS, hardware)?

a workaround is to use arbitrary-precision libraries such as mpmath

It's not a real workaround, because OP want native CPython floats, not mpmath's numbers. In fact, math.ldexp() is used in the mpmath to provide the __float__ dunder for mpf type.

It's possible to use gmpy2:

>>> import gmpy2
>>> gmpy2.set_context(gmpy2.ieee(64))
>>> gmpy2.mpfr(6993274598585239)/gmpy2.mpfr(2)**1126
mpfr('0.0')
>>> gmpy2.set_context(gmpy2.ieee(128))
>>> gmpy2.mpfr(6993274598585239)/gmpy2.mpfr(2)**1126
mpfr('7.67194470417997883016474037150491012e-324',113)
>>> float(_)
1e-323

Though, it's not an option for the mpmath itself.

[tkreuzer]

Just a bug in the platform libc. @tkreuzer, could you please describe your system (OS, hardware)?

Windows 10 on an Intel i/-1065G7.

I implemented a workaround myself now.

def ldexp_manual(x_float, exp):
    """
    Compute x_float * 2**exp for a floating-point number, handling denormals and rounding.

    Args:
        x_float (float): The floating-point number to scale.
        exp (int): The integer exponent of 2 by which to scale.

    Returns:
        float: The result of x_float * 2**exp, respecting IEEE 754 rules.
    """
    # Handle special cases: zero, infinity, or NaN
    if x_float == 0.0 or not math.isfinite(x_float):
        return x_float

    # Get the 64-bit IEEE 754 representation
    bits = struct.unpack('Q', struct.pack('d', x_float))[0]

    # Extract components
    sign = bits >> 63
    exponent = (bits >> 52) & 0x7FF
    mantissa = bits & 0xFFFFFFFFFFFFF

    if exponent == 0:
        # Denormal number
        if mantissa == 0:
            return x_float  # Zero
        # Normalize it
        while mantissa < (1 << 52):
            mantissa <<= 1
            exponent -= 1
        exponent += 1
    else:
        # Normal number, add implicit leading 1
        mantissa |= (1 << 52)

    # Adjust exponent with exp
    new_exponent = exponent + exp

    if new_exponent > 2046:
        # Overflow to infinity
        return math.copysign(math.inf, x_float)
    elif new_exponent <= 0:
        # Denormal or underflow
        if new_exponent < -52:
            # Underflow to zero
            return 0.0

        # Calculate how much we need to shift the mantissa
        mantissa_shift = 1 - new_exponent

        # Get the highest bit, zhat would be shifted off
        round_bit = (mantissa >> (mantissa_shift - 1)) & 1

        # Adjust mantissa for denormal
        mantissa >>= mantissa_shift
        mantissa += round_bit
        new_exponent = 0

    # Reconstruct the float
    bits = (sign << 63) | (new_exponent << 52) | (mantissa & 0xFFFFFFFFFFFFF)
    return float(struct.unpack('d', struct.pack('Q', bits))[0])


def mpf_to_float(mpf_value):
    """
    Convert an mpmath mpf value to the nearest Python float.
    We use ldexp_manual, because ldexp is broken on Windows.
    """
    sign = mpf_value._mpf_[0]
    mantissa = mpf_value._mpf_[1]
    exponent = mpf_value._mpf_[2]

    result = ldexp_manual(mantissa, exponent)
    if sign == 1:
        result = -result
    return result

[skirpichev]

On mpmath's master:

$ python -m timeit -s 'import mpmath;x=mpmath.mpf(1)' 'float(x)'
200000 loops, best of 5: 1.22 usec per loop

With your patch (I doubt it's correct):

$ python -m timeit -s 'import mpmath;x=mpmath.mpf(1)' 'float(x)'
50000 loops, best of 5: 9.09 usec per loop

[tkreuzer]

On mpmath's master:

...

It's slower, I get it. Would not be an issue, if it was implemented in C 🙂

Alternatively you can probably use ldexp for anything that is not a denormal (would need to test that). Then you can simplify the code a bit. And you can replace the shift loop with something like

        # Normalize it
        shift_amount = 52 - mantissa.bit_length()
        mantissa <<= shift_amount
        exponent = 1 - shift_amount

[skirpichev]

Windows 10 on an Intel i/-1065G7.

python -VV output might be helpful.

From quoted issue: Python 3.10.11 (tags/v3.10.11:7d4cc5a, Apr 5 2023, 00:38:17) [MSC v.1929 64 bit (AMD64)] on win32

I would like to see if someone else can reproduce this on same platform.

It's slower, I get it. Would not be an issue, if it was implemented in C

This can be mitigated coding in C, but that's not for the mpmath.

It's also incorrect. Besides raising OverflowError in an appropriate place, you probably should remove if sign == 1 branch in mpf_to_float().

FYI: mpmath/mpmath#945

Alternatively you can probably use ldexp for anything that is not a denormal (would need to test that).

You meant "for anything, that not produces a denormal"? Because 6993274598585239.0 is a normal float. I doubt it's cheap.

CC @tim-one

[tim-one]

The Windows ldexp doesn't round - it chops (truncates: "round to 0"). Whether that's "a bug" is arguable, but Microsoft hasn't budged on this at all.

To make this easier to follow, here's the simplest "failing" case:

>>> import math
>>> tiny = math.nextafter(0, 1)
>>> tiny # smallest non-zero positive float
5e-324
>>> _.hex()
'0x0.0000000000001p-1022'
>>> math.frexp(tiny)
(0.5, -1073) # easier than counting by hand
>>> math.ldexp(0.5, -1073) # ldexp is frexp's inverse
5e-324
>>> math.ldexp(0.75, -1073) # doesn't round up
5e-324
>>> math.ldexp(0.8, -1073) # still doesn't
5e-324
>>> math.ldexp(0.9, -1073) # and still not
5e-324
>>> math.ldexp(0.9999, -1073) # nope!
5e-324
>>> math.ldexp(1.0, -1073) # finally - but not via rounding
1e-323
>>> _.hex()
'0x0.0000000000002p-1022'

So long as the result isn't subnormal, no problem. The original intent of ldexp was to give a fast way to change the exponent without changing the mantissa bits at all. Subnormals don't play nice with that view, though.

Python's math module generally intends to reproduce the platform's C library results for functions of the same name.

There's no clear win to be had here. It's impossible for ldexp to both round and to match C programmers' expectations of how ldexp works on Windows.

If you want the same semantics as multiplying by a power of 2, multiply by a power of 2. The days when ldexp was faster than float multiplication are now ancient history.

Not that it's trivial. In the above, e.g., 2.0**1073 isn't representable as an IEEE double. In bad cases it can take two stages: first use ldexp to shift the input to the smallest normal binade. No information is lost. Then multiply by the appropriate negative power of 2 to finish the job. Let the hardware take care of rounding. There's no need to pick apart - or do arithmetic on - the bits yourself.

[skirpichev]

The Windows ldexp doesn't round - it chops (truncates: "round to 0"). Whether that's "a bug" is arguable

Ah, that explains things for me. Does it happens on all M$ systems or rather limited (e.g. to win32)? (Disclaimer: I had windows box nearby 10+ years ago.)

Sorry, but it seems to be a bug. Standard says: "If a range error due to underflow occurs, the correct result (after rounding) is returned." Rounding mode could be ignored only if no range error occurs.

Python's math module generally intends to reproduce the platform's C library results for functions of the same name.

Are you argue to get rid of CPython's versions of tgamma()/lgamma() and switch to libm's versions? ;-)

So, sometimes we have workarounds for buggy libm functions. This might be the case.

[tim-one]

Does it happens on all M$ systems or rather limited (e.g. to win32)?

I'm not a Windows expert. It's something I've learned from the systems I personally use. It was true under DOS from the start, and is still so under my current 64-bit Windows 10 Pro.

Standard says: "If a range error due to underflow occurs, the correct result (after rounding) is returned." Rounding mode could be ignored only if no range error occur

The great thing about standards is that there are so many of them 😉. You need to identify the exact standard you have in mind. The accepted answer here:

https://stackoverflow.com/questions/32150888/should-ldexp-round-correctly

lets Microsoft off the hook because their C implementation may not define __STDC_IEC_559__ as 1. Well, they didn't know for sure whether it did, but the accepted answer here says __STDC_IEC_559__ is set to 0 on Windows by both gcc and clang (but to 1 on Linux)

https://softwarerecs.stackexchange.com/questions/78793/is-there-any-c-compiler-which-defines-both-stdc-and-stdc-iec-559-to-1

Which isn't the final word either. Hence it's "arguable", but it's not an argument I want to engage in.

Are you argue to get rid of CPython's versions of tgamma()/lgamma() and switch to libm's versions? ;-)

If platform libms have since improved to be better than horrid on those functions, yes, it would be best if we did lose ours. You already recently discovered that, on the platform you used, its gamma is better than ours, We're not doing users on that platform much good by hiding their platform's better implementation. But since we took them over, we have our own backward compatibility worries now.

Microsoft has decades of backward compatibility to worry about with ldexp(), which is why I expect - but don't know - they never change it.

So, sometimes we have workarounds for buggy libm functions. This might be the case.

It might be. But if we do take it over on Windows, our docs need to be clear that Python's ldexp differs in this respect from the native Windows version. The docs should have done so for the gamma functions when we took them over.

It's a fight I don't want to take part in. But I'd be +1 on leaving ldexp() alone, and introducing a new, say, ldexpr() ("ldexp rounded") function that promises to round on all platforms. It would the same as the platform ldexp on most popular platforms.

When I care, I multiply by a negative power of 2. It's cheap to compute 2**-n via ldexp(1.0, -n). But n has to be <= 1074, else that underflows to 0.

>>> math.ldexp(0.75, -1073) # truncates
5e-324
>>> 0.75 * math.ldexp(1, -1073) # rounds
1e-323

So that's the core "trick" to implementing a rounded variant of ldexp on Windows without tedious, slow, and error-prone bit fiddling. frexp() can be used to get the input's exponent. There's no need to look at (let alone change) the mantissa bits.

>>> math.ldexp(6993274598585239, -1126)
5e-324
>>> 6993274598585239 * math.ldexp(1, -1074) * math.ldexp(1., -1126 + 1074)
1e-323

[tim-one]

Here's a very lightly tested workaround for Windows, which defers to the platform function unless rounding may be needed. Rounding (if any is needed) is done by the multiply at the end.

Details

from math import frexp, isfinite, ldexp, nextafter

BITS = 53 # number of mantissa bits
TINY_EXP = frexp(nextafter(0.0, 1.0))[1] # exp of smllest denorm
SMALLEST_NORM_EXP = TINY_EXP + BITS - 1

def ldexpr(x, n):
    if not isfinite(x) or not x:
        # NsNs, infs, and zeros are unchanged.
        return x
    # Non-zero finite.
    if n >= 0:
        # Left shifts may overflow, but don't round.
        return ldexp(x, n)
    # Right shift of non-zero finite.
    original_exp = frexp(x)[1]
    target_exp = original_exp + n
    if target_exp >= SMALLEST_NORM_EXP:
       # No bits are shifted off - no rounding involved.
       return ldexp(x, n)
    # Shifting into denorm-land - trailing bits may be lost.
    if original_exp > SMALLEST_NORM_EXP:
        # Shift down to the smallest normal binade. No bits lost,
        x = ldexp(x, SMALLEST_NORM_EXP - original_exp)
        assert frexp(x)[1] == SMALLEST_NORM_EXP
        n += original_exp - SMALLEST_NORM_EXP
    # Multiplying by 2**n finishes the job, and the HW will round as
    # appropriate, Note; if n < -BITS, all of x is shifted to be <
    # 1/2 ULP of TINY, so should be thrown away. If n is so very
    # negative that ldexp underflows to 0, that's fine; no need to
    # check in advance.
    return x * ldexp(1.0, n)

Note that it's important to do no more than one multiply. lest subtle "double rounding" bugs be introduced. The intent is that no information of any kind is lost before the single multiply.