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

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.

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.

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

[3.15] gh-140036: _pydecimal: avoid slow exponentiation in floor divi…

[3.14] gh-140036: _pydecimal: avoid slow exponentiation in floor divi…

Navigation Menu

Repository navigation

Python implementation of Decimal can hang in floor division #140036

Bug report

Linked PRs

20 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

20 remaining items

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