Skip to content

Commit f3c27a4

Browse files
authored
gh-140036: _pydecimal: avoid slow exponentiation in floor division (#140044)
1 parent 39ea54a commit f3c27a4

3 files changed

Lines changed: 158 additions & 3 deletions

File tree

‎Lib/_pydecimal.py‎

Lines changed: 99 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -447,6 +447,95 @@ def IEEEContext(bits, /):
447447

448448
##### Decimal class #######################################################
449449

450+
# Observation: For all q >= 0 and a >= 1, q < 10**a iff len(str(q)) <= a.
451+
#
452+
# The constants below are used to speed-up "q < 10 ** a" checks to avoid
453+
# computing len(str(q)) as much as possible. Those speed-ups are based on
454+
# the following claims.
455+
#
456+
# See https://github.com/python/cpython/issues/140036 for details.
457+
458+
_LOG_10_BASE_2_LO = float.fromhex('0x1.a934f0979a371p+1')
459+
assert pow(2, _LOG_10_BASE_2_LO) < 10
460+
461+
_LOG_10_BASE_2_HI = float.fromhex('0x1.a934f0979a372p+1')
462+
assert pow(2, _LOG_10_BASE_2_HI) > 10
463+
464+
465+
def _tento(n):
466+
"""Compute 10 ** n with 1 base-5 exponentiation and 1 bit-shift."""
467+
return (5 ** n) << n
468+
469+
470+
def _is_less_than_pow10a_use_str(q, a):
471+
"""Try to efficiently check len(str(q)) <= a, or equivalently q < 10**a.
472+
473+
If the comparison cannot be obtained from q.bit_length(),
474+
then str(q) is explicitly computed and may raise ValueError.
475+
476+
Return (len(str(q)) <= a, None) or (len(str(q)) <= a, str(q)).
477+
"""
478+
if q.bit_length() < a * _LOG_10_BASE_2_LO:
479+
# Claim: If 0 < z <= log2(10) and q.bit_length() < a*z, then q < 10**a.
480+
# Proof: By contradiction, q >= 10**a. By definition,
481+
# log2(q) >= a*log2(10) >= a*z > q.bit_length().
482+
# In particular, q > 2**q.bit_length(), which is impossible.
483+
484+
# assert q < 10 ** a
485+
return True, None
486+
elif q.bit_length() >= 1 + a * _LOG_10_BASE_2_HI:
487+
# Claim: If z > log2(10) and q.bit_length() >= 1 + a*z, then q > 10**a.
488+
# Proof: Since q >= 2**(q.bit_length()-1), we have
489+
# q >= 2**(q.bit_length()-1) >= 2**(a*z) > 2**(a*log2(10)) = 10**a.
490+
491+
# assert q > 10 ** a
492+
return False, None
493+
# Handle cases that fail due to floating point precision loss
494+
# when computing _LOG_10_BASE_2_LO and _LOG_10_BASE_2_HI, or
495+
# that cannot be distinguished with (q.bit_length(), a) only.
496+
#
497+
# For instance, (q1, a) = (95, 2) and (q2, a) = (105, 2) produce
498+
# different results but q1.bit_length() == q2.bit_length() == 7.
499+
str_q = str(q) # can raise a ValueError
500+
is_valid = len(str_q) <= a
501+
return is_valid, str_q
502+
503+
504+
def _is_less_than_pow10a(q, a, *, exact=True, ulp_order=20):
505+
"""Check that len(str(q)) <= a without computing str(q).
506+
507+
When *exact* is false, computing len(str(q)) is replaced by f(q):
508+
509+
f(q) = floor(log10(q) + ulp(log10(q)) * ulp_order + 1.0)
510+
511+
Most of the time, f(q) = len(str(q)) but in some cases, it may
512+
happen that f(q) > len(str(q)).
513+
514+
When *exact* is true, computing len(str(q)) requires one bigint
515+
exponentiation that only depends on q.
516+
"""
517+
if q < 10:
518+
return a >= 1
519+
520+
z = _math.log10(q)
521+
t = _math.ulp(z) * ulp_order
522+
523+
if exact:
524+
intlo = int(z - t)
525+
inthi = int(z + t)
526+
diff = inthi - intlo
527+
assert diff in (0, 1)
528+
if diff == 1:
529+
lo = _tento(inthi) # may be slow
530+
if q < lo:
531+
inthi -= 1
532+
assert q >= (lo // 10)
533+
ndigits = inthi + 1
534+
else:
535+
ndigits = int(z + t + 1.0)
536+
return ndigits <= a
537+
538+
450539
# Do not subclass Decimal from numbers.Real and do not register it as such
451540
# (because Decimals are not interoperable with floats). See the notes in
452541
# numbers.py for more detail.
@@ -1360,8 +1449,13 @@ def _divide(self, other, context):
13601449
else:
13611450
op2.int *= 10**(op2.exp - op1.exp)
13621451
q, r = divmod(op1.int, op2.int)
1363-
if q < 10**context.prec:
1364-
return (_dec_from_triple(sign, str(q), 0),
1452+
is_valid, str_q = _is_less_than_pow10a_use_str(q, context.prec)
1453+
if is_valid:
1454+
if str_q is None:
1455+
str_q = str(q)
1456+
# assert q < 10 ** context.prec
1457+
# assert len(str(q)) <= context.prec
1458+
return (_dec_from_triple(sign, str_q, 0),
13651459
_dec_from_triple(self._sign, str(r), ideal_exp))
13661460

13671461
# Here the quotient is too large to be representable
@@ -1516,7 +1610,9 @@ def remainder_near(self, other, context=None):
15161610
r -= op2.int
15171611
q += 1
15181612

1519-
if q >= 10**context.prec:
1613+
if not _is_less_than_pow10a(q, context.prec):
1614+
# assert q >= 10 ** context.prec
1615+
# assert len(str(q)) > context.prec
15201616
return context._raise_error(DivisionImpossible)
15211617

15221618
# result has same sign as self unless r is negative

‎Lib/test/test_decimal.py‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@
2424
with the corresponding argument.
2525
"""
2626

27+
import contextlib
2728
import logging
2829
import math
2930
import os, sys
@@ -2612,6 +2613,37 @@ def tearDown(self):
26122613
sys.set_int_max_str_digits(self._previous_int_limit)
26132614
super().tearDown()
26142615

2616+
def test_helper__is_less_than_pow10a_use_str_slow_path(self):
2617+
# Test the "slow" path of _is_less_than_pow10a_use_str().
2618+
a, b = 2, 7
2619+
2620+
# Choose q1, q2 such that len(str(q1)) <= a < len(str(q2))
2621+
# and q1.bit_length() == q2.bit_length() == b to check that
2622+
# we cover the "slow" path correctly even for small values.
2623+
q1, q2 = 95, 105
2624+
b1, b2 = q1.bit_length(), q2.bit_length()
2625+
2626+
self.assertEqual(b1, b)
2627+
self.assertEqual(b2, b)
2628+
2629+
# ensure that the first "fast" check doesn't hold
2630+
self.assertGreaterEqual(b, a * self.decimal._LOG_10_BASE_2_LO)
2631+
# ensure that the second "fast" check doesn't hold
2632+
self.assertLess(b, 1 + a * self.decimal._LOG_10_BASE_2_HI)
2633+
2634+
cond_q1, str_q1 = self.decimal._is_less_than_pow10a_use_str(q1, a)
2635+
self.assertTrue(cond_q1)
2636+
self.assertIsNotNone(str_q1)
2637+
2638+
cond_q2, str_q2 = self.decimal._is_less_than_pow10a_use_str(q2, a)
2639+
self.assertFalse(cond_q2)
2640+
self.assertIsNotNone(str_q2)
2641+
2642+
def test_helper__is_less_than_pow10a(self):
2643+
# TODO(picnixz): find a simple test case with custom ulp_order.
2644+
pass
2645+
2646+
26152647
class PythonAPItests:
26162648

26172649
def test_abc(self):
@@ -4556,6 +4588,15 @@ def test_decimal_attributes(self):
45564588

45574589
class Coverage:
45584590

4591+
@contextlib.contextmanager
4592+
def unbound_context(self, prec=None, Emax=None, Emin=None):
4593+
with self.decimal.localcontext() as c:
4594+
c.prec = self.decimal.MAX_PREC if prec is None else prec
4595+
c.Emax = self.decimal.MAX_EMAX if Emax is None else Emax
4596+
c.Emin = self.decimal.MIN_EMIN if Emin is None else Emin
4597+
c.traps[self.decimal.Inexact] = 1
4598+
yield c
4599+
45594600
def test_adjusted(self):
45604601
Decimal = self.decimal.Decimal
45614602

@@ -4723,6 +4764,22 @@ def test_divmod(self):
47234764
self.assertTrue(c.flags[InvalidOperation] and
47244765
c.flags[DivisionByZero])
47254766

4767+
def test_divide_unbound_context(self):
4768+
with self.unbound_context() as c:
4769+
x = self.decimal.Decimal('1')
4770+
y = x // 1 # should be fast
4771+
4772+
def test_remainder_near(self):
4773+
L = 1000
4774+
limit = sys.get_int_max_str_digits()
4775+
sys.set_int_max_str_digits(L)
4776+
self.addCleanup(sys.set_int_max_str_digits, limit)
4777+
4778+
with self.unbound_context(prec=2 * L) as c:
4779+
self.assertEqual(c.prec, 2 * L)
4780+
x = self.decimal.Decimal(f'1e{L}')
4781+
y = x.remainder_near(1) # must not raise a ValueError
4782+
47264783
def test_power(self):
47274784
Decimal = self.decimal.Decimal
47284785
localcontext = self.decimal.localcontext
Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
Avoid hanging in floor division of pure Python :class:`decimal.Decimal`
2+
instances when the context precision is very large. Patch by Bénédikt Tran.

0 commit comments

Comments
 (0)