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

I cannot say why but I fear that we could have corner cases where the approximation would not be enough. An alternative I had in mind was to check len(str(q)) <= context.prec. Note that, in _divide we already call str(q) if the precision check passes. If str(q) fails afterwards, we'll also have a ValueError. So we need str(q) to succeed to get a meaningful division.

In particular, I don't know how we should handle errors with str(q) when q is very large. So I think we can just use str(q) (there are many other places where str(...) is used without bothering about the int2digits limit).

EDIT: we can just make a first rough check (yours) and then recheck with len(str(...)) as we anyway need to compute str(q)!

An alternative I had in mind was to check len(str(q)) < context.prec.

You want <= rather than <, but I agree this is the way to go (with some refactoring to avoid computing str(q) twice on the non-error path).

There's similar code in remainder_near that should be fixed, too:

cpython/Lib/_pydecimal.py

Lines 1514 to 1515 in 52996aa

    
           if q >= 10**context.prec: 
        
               return context._raise_error(DivisionImpossible)

.

I'm a bit worried here because, unlike in the previous case, we don't need to compute str(q) in the success path. For instance:

sys.set_int_max_str_digits(640)
...
_unbounded_dec_context.prec = 1000
x = decimal.Decimal('1e640')
y = x.remainder_near(1)

With the main version, this is fast. But with the patched version relying on str we would have a ValueError because of the conversion. So I don't know how we should handle this case (for now, I'm going to only fix the first case and not this one)

But with the patched version relying on str we would have a ValueError because of the conversion.

Yep, I think there's a ton of int-to-str conversion issues lurking here.

Yep, I think there's a ton of int-to-str conversion issues lurking here.

Python 3.13.7 (main, Aug 15 2025, 23:01:14) [Clang 17.0.0 (clang-1700.0.13.5)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from _pydecimal import Decimal, getcontext
>>> getcontext().prec = 10000
>>> Decimal(2) / Decimal(3)
Traceback (most recent call last):
  File "<python-input-2>", line 1, in <module>
    Decimal(2) / Decimal(3)
    ~~~~~~~~~~~^~~~~~~~~~~~
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.13/lib/python3.13/_pydecimal.py", line 1293, in __truediv__
    ans = _dec_from_triple(sign, str(coeff), exp)
                                 ~~~^^^^^^^
ValueError: Exceeds the limit (4300 digits) for integer string conversion; use sys.set_int_max_str_digits() to increase the limit
>>> Decimal(2).ln()
Traceback (most recent call last):
  File "<python-input-4>", line 1, in <module>
    Decimal(2).ln()
    ~~~~~~~~~~~~~^^
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.13/lib/python3.13/_pydecimal.py", line 3153, in ln
    coeff = _dlog(c, e, places)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.13/lib/python3.13/_pydecimal.py", line 5758, in _dlog
    log_d = _ilog(c, 10**p) # error < 5 + 22 = 27
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.13/lib/python3.13/_pydecimal.py", line 5692, in _ilog
    T = -int(-10*len(str(M))//(3*L))
                     ~~~^^^
ValueError: Exceeds the limit (4300 digits) for integer string conversion; use sys.set_int_max_str_digits() to increase the limit
>>> Decimal(2).sqrt()
Traceback (most recent call last):
  File "<python-input-5>", line 1, in <module>
    Decimal(2).sqrt()
    ~~~~~~~~~~~~~~~^^
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.13/lib/python3.13/_pydecimal.py", line 2730, in sqrt
    ans = _dec_from_triple(0, str(n), e)
                              ~~~^^^
ValueError: Exceeds the limit (4300 digits) for integer string conversion; use sys.set_int_max_str_digits() to increase the limit

What should we do here? should we temporarily raise the limit for those conversions? (it would be funny for people who want infinite precision/exact computations to be actually limited by the conversion limit...).

What should we do here?

Very good question. :-)

Should we temporarily raise the limit for those conversions?

I don't think that's an option while all we have is a global setting in sys, since changing the limit in one thread changes it in all threads (so a thread running Decimal operations could inadvertently cause changed behaviour in some other thread).

Ah. I thought it was a per-thread limit not a global one. Can't we just expose on int the number of digits that we get through _PyLong_DigitCount? Alternatively, we could expose it as a private function in _math just to avoid importing other modules.

the number of digits that we get through _PyLong_DigitCount?

That's not what you want, I'm afraid: that's the number of internal base-power-of-two digits (base 2**15 or 2**30 depending on configure-time settings). I don't think there's anything in the codebase that reliably computes the exact number of decimal digits of a Python int without actually computing the str.

There's similar code in remainder_near that should be fixed

There might be more. Some code use exponents in same way, which might be enormous as well (MAX_EMAX/MIN_EMIN).

Yep, I think there's a ton of int-to-str conversion issues lurking here.

And I think this deserve own issue.

What should we do here? should we temporarily raise the limit for those conversions?

IMO, in a long term - use better bigint library and avoid the need of this limit ;-)

This limitation should be set by application. I.e. no limits by default. And, maybe, API could be better with a context manager instead of some global setting.

Can't we just expose on int the number of digits that we get through _PyLong_DigitCount?

We have bit_length() method. Assuming this and known layout of the int (available in sys) - it's easy to calculate number of "digits". But why?

We have bit_length() method. Assuming this and known layout of the int (available in sys) - it's easy to calculate number of "digits". But why?

That was for the check q > 10**prec where we don't want to replace by len(str(q)). But for other places, it wouldn't be helpful. (So what I wanted here is just a very fast len(str(q)) that isn't limited by sys.get_int_max_str_digits()).

@mdickinson

@mdickinson

I don't think there's anything in the codebase that reliably computes the exact number of decimal digits of a Python int without actually computing the str

Maybe we should address that directly. Here's a simple approach that arrives at the result iteratively, relying on that $\log_2{10}$ is bounded above by 4.

https://stackoverflow.com/questions/18828377/biginteger-count-the-number-of-decimal-digits-in-a-scalable-method/18860385#18860385

The loop can require very large ints indeed,, but they shrink quickly across iterations, and even the first iteration may be able to exploit CPython's newer sub-quadratic bigint division.

I don't believe a "truly fast" way exists that's always exact. Proof: if one did, GMP wouldn't have given up on it 😉.

BTW, I haven't implemented or tested that approach. It just "feels like it's probably on track" to me.

The loop can require very large ints indeed, but they shrink quickly across iterations

My quick tests:

Benchmark	ndigits	log10
10**1	573 ns	373 ns: 1.54x faster
10**10	685 ns	1.34 us: 1.95x slower
10**100	1.93 us	4.58 us: 2.37x slower
10**1000	42.5 us	25.6 us: 1.66x faster
10**10000	1.29 ms	879 us: 1.47x faster
10**100000	43.5 ms	30.1 ms: 1.45x faster
Geometric mean	(ref)	1.03x faster

Details

# bench.py

import os
import sys

sys.set_int_max_str_digits(0)

import pyperf

def ndigits(n):
    return len(str(n))

def log10(n):
    ndigits = 1
    while (nbits := n.bit_length()) > 4:
        reduce = nbits // 4;
        n //= 10**reduce
        ndigits += reduce
    return ndigits + (n > 9)

f = eval(os.getenv('_F'))

runner = pyperf.Runner()
for i in range(6):
    e = 10**i
    n = 10**e
    h = f'10**{e}'
    runner.bench_func(h, f, n)

Proof: if one did, GMP wouldn't have given up on it

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.

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.

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…

gh-140036: Test modulo and divmod of pure Python Decimal in unbounded…

I've asked Claude to investigate how to fix the modulo operation but there were other cases that could likely be slow. More generally, I want to see if we can defer 10**k computations when k is very large as much as possible.

[3.15] gh-140036: Test modulo and divmod of pure Python Decimal in un…

[3.14] gh-140036: Test modulo and divmod of pure Python Decimal in un…

Navigation Menu

Repository navigation

Python implementation of Decimal can hang in floor division #140036

Bug report

Linked PRs

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

	if q >= 10**context.prec:
	return context._raise_error(DivisionImpossible)

Python implementation of Decimal can hang in floor division #140036

Description

Bug report

Linked PRs

Activity

picnixz commented on Oct 13, 2025

mdickinson commented on Oct 13, 2025

mdickinson commented on Oct 13, 2025

picnixz commented on Oct 13, 2025

mdickinson commented on Oct 13, 2025

mdickinson commented on Oct 13, 2025

picnixz commented on Oct 13, 2025

mdickinson commented on Oct 13, 2025

picnixz commented on Oct 13, 2025

mdickinson commented on Oct 13, 2025

skirpichev commented on Oct 13, 2025

picnixz commented on Oct 13, 2025

tim-one commented on Oct 14, 2025

skirpichev commented on Oct 14, 2025

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

picnixz commented on Oct 28, 2026

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions