Python implementation of Decimal can hang in floor division

@tim-one

Bug report

# import decimal
import _pydecimal as decimal

_unbounded_dec_context = decimal.getcontext().copy()
_unbounded_dec_context.prec = decimal.MAX_PREC
_unbounded_dec_context.Emax = decimal.MAX_EMAX
_unbounded_dec_context.Emin = decimal.MIN_EMIN
_unbounded_dec_context.traps[decimal.Inexact] = 1 # sanity check
decimal.setcontext(_unbounded_dec_context)

x = decimal.Decimal('1')
y = x // 1

The C implementation evaluates the last line instantly, but the Python implementation hangs for a long time (infinitely?). The modulus operation has the same behavior.

cc @tim-one

Linked PRs

gh-140036: _pydecimal: avoid slow exponentiation in floor division #140044

In _divide, it tries to compute 10 ** context.prec which in this case is 10 ** 999999999999999999 to compare to q. However, taking the base-10 logarithm on both sides may result in a precision loss and have an incorrect result.

In _divide, it tries to compute 10 ** context.prec which in this case is 10 ** 999999999999999999 to compare to q. However, taking the base-10 logarithm on both sides may result in a precision loss and have an incorrect result.

But why it's possible to calculate this in the C implementation?
Those implementations aren't the same in terms of precision?

It's possible because the backend actually stores the number of digits of the divmod result. So you can just do some q->digits > ctx->prec check without having to re-raise to the power of 10. (I believe the implementations are identical for the precision handling; it just that the C implementation doesn't have some $10^{prec}$ check because it knows the number of digits in q already while the Python implementation doesn't).

Maybe just:

if q <= 0 or _math.log10(q) < context.prec:
    ...

We could speed-up some cases by relying on .bit_length() to get the maximum number of the next power of 2. We can switch to that estimation if we know in advance that 10**context.prec will take a long time (for instance, on my machine 10 ** 1000000 becomes slow, so anything below can be directly computed).

if q <= 0 or _math.log10(q) < context.prec:

For most cases, it would be fine, but _math.log10 may either overestimate or underestimate if we exceed 15 digits. So I don't think it would sufficient here.

Oh, a similar issue was already discussed in #118027 and #118164.

For most cases, it would be fine, but _math.log10 may either overestimate or underestimate if we exceed 15 digits.

I assume special handling of integers, available in math.log*():

cpython/Modules/mathmodule.c

Lines 2302 to 2313 in 52996aa

    
                   if (x == -1.0 && PyErr_Occurred()) { 
        
                       if (!PyErr_ExceptionMatches(PyExc_OverflowError)) 
        
                           return NULL; 
        
                       /* Here the conversion to double overflowed, but it's possible 
        
                          to compute the log anyway.  Clear the exception and continue. */ 
        
                       PyErr_Clear(); 
        
                       x = _PyLong_Frexp((PyLongObject *)arg, &e); 
        
                       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; 
        
                   }

Wouldn't we still have a problem with the exact computation of the log as the result is a float?
For instance, with a context precision of 30 and x = decimal.Decimal('1e30') - 1, we have $q = 999999999999999999999999999999$. However math.log10(q) == 30 so the condition log10(q) < ctx.prec would fail, whereas q < 10 ** 30 would succeed (that's where we underestimate the log10 because of the inexact return type)

You are right, it should be something like q.bit_length() < 1 + context.prec*3.32192809488736.

Ya, that's the idea. For "small enough" arguments, probably better to do len(str(n)) directly rather than enter the loop.

I doubt it's useful enough to warrant some API in GMP. It already has mpn_sizeinbase(), which result is either exact or 1 too big.

Please. This is GMP 😉. There is no heroic effort too daunting for them to undertake, no matter how useless. Instead their docs say that if xxx_sizeinbase() isn't good enough, convert it to a string and get its length. If they knew of any "reallly fast" way to get it always right, they almost certainly would have implemented it, regardless of implementation effort.

So the other question here is whether "always right" is actually needed. Is "exact or 1 too big" actually good enough? It is, e.g., if you only need to allocate enough space to hold the string. "One too big" is so insignificant in that context it's not worth any cycles to avoid.

Instead their docs say that if xxx_sizeinbase() isn't good enough, convert it to a string and get its length.

Hmm, that's something I miss for mpz_sizeinbase()/mpn_sizeinbase(). (Just did a quick search.) Maybe it was present in older manuals, I'm looking on v6.3.0.

if you only need to allocate enough space to hold the string.

IMO, it's most common task for this function for base=10, which hardly warrant more implementation efforts (and computational cost) to be exact.

Hmm, that's something I miss for mpz_sizeinbase()/mpn_sizeinbase(). (Just did a quick search.) Maybe it was present in older manuals, I'm looking on v6.3.0.

I may well have conflated that with Paul Zimmermann's reply here, which is certainly part of the "unpublished docs" 😉:

https://gmplib.org/list-archives/gmp-bugs/2020-August/004859.html

If you want the exact size, use mpz_get_str() and strlen(), or mpz_ui_pow_ui() and mpz_cmp().

For the rest, I'm not arguing about it. I think it dead obvious that if they knew of a very fast, always right way, they would have implemented it. You don't. In context, "the truth" makes no difference to this issue report.

If "exact or 1 too big" is good enough for our purposes (which I don't know), then that's something we could pursue to get an answer much faster.

If we need exact answers, try something like this, which does at worst only one "big power". Fast for small inputs, slower as they get larger.

from math import log10
def ndig(n):
    assert n >= 0
    if n < 10:
        return 1
    trial = int(log10(n)) + 1 # always fast
    lo = 10 ** (trial - 1)
    if n < lo:
        trial -= 1
        lo //= 10
        assert n >= lo
    else:
        hi = lo * 10
        if n >= hi:
             trial += 1
             assert n < hi * 10
     return trial

The asserts take the place of rigorous error analysis 😉. But I don't want to take the time for that unless there's real interest in pursuing this. The heart of it is simply that log10() "is close", pretty much regardless of how sloppy (within reason) the platform libm is, and that Python's implementation of log10(bigint) only looks at the most-significant ~53 bits of bigint.

And another. No wrong results seen so far in testing, but don't yet have high confidence. "Features":

Usually gets away with no bigint ops apart from 1 shift and 1 add.
In "close to a power of 10" cases, resolved with a single bigint pow,
But the tail including that part can be commented out if you're willing to settle for "exact or 1 larger"
Get rid of the print() statements 😉.

from math import log10, ulp

def ndig(n):
    assert n >= 0
    if n < 10:
        return 1
    # Get an upper bound.
    m = n + (n >> 52) # by boosting the input
    raw = log10(m)
    # And by boosting the computed log. Hard to believe a modern
    # libm would be off by 20 ULP,
    # Note that both of these boosts are very small in relative terms.
    raw += ulp(raw) * 20
    raw += 1.0
    trial = int(raw)
    # The rest can be skipped if willing to settle for maybe one too
    # large.
    frac = raw - trial
    if frac < 0.04:  # a little above a power of 10
        print("?", end="")
        lo = 10 ** (trial - 1)
        if n < lo:
            print("-", end="")
            trial -= 1
            lo //= 10
            assert n >= lo
    else:
        print('=', end="")
    return trial

Usually gets away with no bigint ops apart from 1 shift and 1 add.

Which appear unneeded. Those operations rarely changed the computed log, and seemingly never by more than 1 ULP (which I believe is provable - log10 grows too slowly for such a small delta to matter much).

So

m = n + (n >> 52) # by boosting the input
raw = log10(m)

can be replaced by plain raw = log10(n) without harm.

Other tweaks:

10**n is faster spelled as (5**n) << n.
0.04 could certainly be made smaller. How small? Still unsure.
Error analysis looks easy enough to finish if there are no more than a billion decimal digits. At a trillion, though, we're consuming about 40 bits of the 53 we get back just to pick out the integer part. It gets dicier then. If that ever becomes "a problem", it would be straightforward (just annoying) to use decimal itself to compute the log to higher precision. A bonus is that decimal.log10() promises to return a correctly rounded result (while libm promises nothing about max error).

@picnixz

@picnixz's proposal was for int-level API for this. I doubt it does make sense. Looks too low-level, while we don't have even a method to print integer in a generic base (supported by the int's constructor).

But the _pydecimal is a very specific beast. Here it seems to be a common story to query len(str(someint)) (grep shows ~40 cases vs ~80 str() calls). Few others use cases might come from comparisons like in the current issue (someint op 10**hugepower, where later come from big precision or exponent).

BTW, the libtommath has exact mpz_sizeinbase() analog (called mp_radix_size(), based on slow integer logarithm). Though, their mp_radix_size_overestimate() have much bigger error bound (22!) than the GMP (1).

Yes my proposal was to just have a function that computes len(str(n)) as fast as possible, without necessarily having to compute str(n) itself just for the decimal module.

Is it OK if it's one too large at times? I gave a function above that takes small time independent of n's size if that's OK. But if exact is needed, it can (but usually doesn't) need a possibly expensive bigint pow() to subtract 1 when needed.

It depends how often it will be wrong (I had the same issue for reprlib...). If we overestimate the length, we will only say "we cannot do the division because the precision is too low" even if it is not the case. If this is acceptable, then maybe.

For now I only focused on avoiding n < 10^p as much as possible and we were able to find a nice formula for that which should suffice most of the time (though we may want to compute the number of cases it succeeds to give the answer without computing len(str(n))). However, as Sergey observed, we have other cases of len(str(n)) and some of them may require the exact number of digits.

Whatever you do, I strongly advise creating a distinct function to hold the logic, rather than duplicate code inline. Things can be surprising. For example, you're implementing a form of interval arithmetic by using two constants to bracket the true value of $\log_{2}{10}$ But the multiplies use nearest/even rounding, not the $-\infty$ and $+\infty$ rounding appropriate for actual interval arithmetic. Maybe that shortcut matters in some cases, maybe not. Or, for another, context precision can be set to a value that requires more than 53 bits to express, in which case any floating-point arithmetic using it is blind to the lowest-order bits.

Rigorous error analysis is rarely easy. My own function becomes dicey when the number of decimal digits heads into the trillions, for reasons having to do with that 40 bits eats up so much of the computed value just to get the integer part that the room to absorb floating-point "slop" gets uncomfortably low. I could increase the headroom by adopting (as you did) hard-coded constants rather than cater to worst-case platform log10() implementations, at some cost in additional code complexity.

But a trillion decimal digits is unrealistic for current HW and plausible apps. By putting it all in a single function, if/when that changes there's an obvious, single place to make changes.

we may want to compute the number of cases it succeeds to give the answer without computing len(str(n)

The function I gave never converts to string, but gets the same effect (when needed) by comparing against a bigint 10**k (generally faster than converting to string if there are lots of digits)

Assuming all ints in range(10**k, 10**(k+1)) are equally likely, the given 0..04 ends up needing to check about 1 in 90 cases ($0,04 \approx \log_{10}{1.1}$). Changing that to 0.004 ($\approx \log_{10}{1.01}$) would reduce that to about 1 in 900, but puts more pressure on getting a more accurate log10. What's there seems fine for "merely" millions of decimal digits, but may get in trouble with billions (haven't tried to analyze that sharply enough, but testing hasn't found a problem yet).

Of course it's not possible to run exhaustive tests on cases so large; even checking just hundreds can take hours of CPU time.

And one more, using a crude but effective way of using interval arithmetic. No check is needed unless an integer is in the interior of a log10 interval. This is sensitive enough so that, e.g., passing 10**9 does need to check, but passing one more or less than that doesn't. An input has to be "very close" to a power of 10 to trigger a check.

math.log10(bigint) is always fast. It effectively uses just the leading 53 bits of the input, plus its power-of-2 exponent. It does more than actually needed here (the leading 53 bits are subject to nearest/even rounding),.
It detects on its own when machine precision has reached its limits. As is, that requires an input with over 250 trillion decimal digits, so not a realistic concern for now.
It would be easy to add on optional exact=True argument. When False, no bigint pow would ever be needed.
The print() and tags[] things are there to help testing, and serve no algorithmic purpose.
$\pm 20$ ULP is intended to be conservative, but you never now how bad a libm log10 could be.

from math import log10, ulp

def tento(n):
    return (5 ** n) << n

tags = [0] * 3

def ndig(n):
    assert n >= 0
    if n < 10:
        tags[0] += 1
        return 1
    raw = log10(n) # always fast
    delta = ulp(raw) * 20
    intlo = int(raw - delta)
    inthi = int(raw + delta)
    diff = inthi - intlo
    assert diff >= 0
    tag = 0
    #print('diff', diff, raw, delta, intlo, inthi)
    if not diff:
        print('=', end="")
    elif diff == 1:
        tag = 1
        print("?", end="")
        lo = tento(inthi)
        if n < lo:
            print("-", end="")
            tag = 2
            inthi -= 1
            lo //= 10
            assert n >= lo
    else:
        raise SystemError(diff) # flesh this out
    tags[tag] += 1
    return inthi + 1

But the multiplies use nearest/even rounding, not the − ∞ and + ∞ rounding appropriate for actual interval arithmetic

Oh great. That's why I always hated rounding issues because of all the different possibilities!

Or, for another, context precision can be set to a value that requires more than 53 bits to express

Urgh, yes, the max prec is 999999999999999999, which is roughly $2^{60}$.

Whatever you do, I strongly advise creating a distinct function to hold the logic, rather than duplicate code inline

I think I'll have two functions, one that will eventually compute len(str(n)) as it's needed in later computations anyway and one that doesn't need to compute it using the ndig function you suggested (I also had a local version of that idea, though the check was a bit different). The problem I had with this is that it entirely depended on log10 implementation which I wanted to avoid.

Now, let's be realisitic. First of all, this only happens if someone is using the pure Python implementation. Most of the time, I think our off-by-one mark could be fine. The check is $len(str(n)) \leq precision \Longrightarrow OK$. We have two cases:

Case 1: str(n) must be explicitly computed if the check passes. No need for a contrived function, we can just compute str(n) and take len(str(n)). Here, we will never have an issue.
Case 2: str(n) does not need to be computed. We can overestimate len(str(n)) by some ndig(n) assuming that $ndig(n) \geq len(str(n))$ with equality most of the time.

In case 2, overestimating means that we're actually performing a stricter check so it's fine. It would be problematic if we underestimated the number of digits. So I think I'm just going with your initial suggestion and leave it with overestimation if it's about bound checking to see if a computation works or not.

Navigation Menu

Repository navigation

Python implementation of Decimal can hang in floor division #140036

Bug report

Linked PRs

14 remaining items

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

	if (x == -1.0 && PyErr_Occurred()) {
	if (!PyErr_ExceptionMatches(PyExc_OverflowError))
	return NULL;
	/* Here the conversion to double overflowed, but it's possible
	to compute the log anyway. Clear the exception and continue. */
	PyErr_Clear();
	x = _PyLong_Frexp((PyLongObject *)arg, &e);
	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;
	}

Python implementation of Decimal can hang in floor division #140036

Description

Bug report

Linked PRs

Activity

picnixz commented on Oct 13, 2025

efimov-mikhail commented on Oct 13, 2025

picnixz commented on Oct 13, 2025

skirpichev commented on Oct 13, 2025

picnixz commented on Oct 13, 2025

serhiy-storchaka commented on Oct 13, 2025

skirpichev commented on Oct 13, 2025

picnixz commented on Oct 13, 2025

skirpichev commented on Oct 13, 2025

14 remaining items

tim-one commented on Oct 14, 2025

skirpichev commented on Oct 14, 2025

tim-one commented on Oct 14, 2025

tim-one commented on Oct 14, 2025

tim-one commented on Oct 14, 2025

tim-one commented on Oct 15, 2025

tim-one commented on Oct 15, 2025

skirpichev commented on Oct 15, 2025

picnixz commented on Oct 16, 2025

tim-one commented on Oct 16, 2025

picnixz commented on Oct 16, 2025

tim-one commented on Oct 16, 2025

tim-one commented on Oct 17, 2025

tim-one commented on Oct 20, 2025

picnixz commented on Oct 26, 2025

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions