Python implementation of Decimal can hang in floor division

[serhiy-storchaka]

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

[picnixz]

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.

[efimov-mikhail]

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?

[picnixz]

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).

[skirpichev]

Maybe just:

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

[picnixz]

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.

[serhiy-storchaka]

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

[skirpichev]

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;
	}

[picnixz]

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)

[skirpichev]

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