gh-118750: Asymptotically faster int(string)

Lib/_pylong.py

@@ -210,6 +210,121 @@ def inner(a, b):
     w5pow = compute_powers(len(s), 5, DIGLIM)
     return inner(0, len(s))
 
+# Asymptotically faster version, using the C decimal module. See
+# comments at the end of the file. This uses decimal arithmetic to
+# convert from base 10 to base 256. The latter is just a string of
+# bytes, which CPython can convert very efficiently to a Python int.
+
+# log of 10 to base 256 with best-possible 53-bit precision. Obtained
+# via:
+#    from mpmath import mp
+#    mp.prec = 1000
+#    print(float(mp.log(10, 256)).hex())
+_LOG_10_BASE_256 = float.fromhex('0x1.a934f0979a371p-2') # about 0.415
+
+def _dec_str_to_int_inner(s):
+    BYTELIM = 512
+    D = decimal.Decimal
+    result = bytearray()
+
+    def inner(n, w):
+        #assert n < D256 ** w # required, but too expensive to check
+        if w <= BYTELIM:
+            # XXX Stefan Pochmann discovered that, for 1024-bit ints,
+            # `int(Decimal)` took 2.5x longer than `int(str(Decimal))`.
+            # So simplify this code to the former if/when that gets
+            # repaired.
+            result.extend(int(str(n)).to_bytes(w)) # big-endian default
+            return
+        w2 = w >> 1
+        if 0:
+            # This is maximally clear, but "too slow". `decimal`
+            # division is asymptotically fast, but we have no way to
+            # tell it to reuse the high-precision reciprocal it computes
+            # for pow256[w2], so it has to recompute it over & over &
+            # over again :-(
+            hi, lo = divmod(n, pow256[w2][0])
+        else:
+            p256, recip = pow256[w2]
+            # The integer part will have a number of digits about equal
+            # to the difference between the log10s of `n` and `pow256`
+            # (which, since these are integers, is roughly approximated
+            # by `.adjusted()`). That's the working precision we need,
+            # but add some guard digits to protect against the "about"
+            # and "roughly" uncertainties.
+            ctx.prec = max(n.adjusted() - p256.adjusted(), 0) + 8
+            hi = +n * +recip # unary `+` chops back to ctx.prec digits
+            ctx.prec = decimal.MAX_PREC
+            hi = hi.to_integral_value() # lose the fractional digits
+            lo = n - hi * p256
+            # Because we've been uniformly rounding down, `hi` is a
+            # lower bound on the correct quotient.
+            assert lo >= 0
+            # Adjust quotient up if needed. It usually isn't. In random
+            # testing, the loop body entered about one in 100 thousand
+            # cases. I never saw it need more than one iteration.
+            count = 0
+            while lo >= p256:
+                count += 1
+                # If the assert fails, chances are decent we're sooooo
+                # far off it may seem to run forever otherwise - the
+                # error analysis was fatally flawed in this case.
+                assert count < 10, (count, w, len(s),
+                                    n.adjusted(), p256.adjusted())
+                lo -= p256
+                hi += 1
+            # The assert should always succeed, but way too slow to keep
+            # enabled.
+            #assert hi, lo == divmod(n, pow256[w2][0])
+        inner(hi, w - w2)
+        inner(lo, w2)
+
+    # How many base 256 digits are needed?. Mathematically, exactly
+    # floor(log256(int(s))) + 1. There is no cheap way to compute this.
+    # But we can get an upper buond, and that's necessary for our error
+    # analysis to make sense. int(s) < 10**len(s), so the log needed is
+    # < log256(10**len(s)) = len(s) * log256(10). However, using
+    # finite-precision floating point for this, it's possible that the
+    # computed value is a little less than the true value. If the true
+    # value is at - or a little higher than - an integer, we can get an
+    # off-by-1 error too low. So we add 2 instead of 1 if chopping lost
+    # a fraction > 0.9.
+    log_ub = len(s) * _LOG_10_BASE_256
+    log_ub_as_int = int(log_ub)
+    w = log_ub_as_int + 1 + (log_ub - log_ub_as_int > 0.9)
+    # And what if we'vv plain exhausted the limits of HW floats? We
+    # could compute the log to any desired precision using `decimal`,
+    # but it's not plausible that anyone will pass a string requiring
+    # trillions of bytes (unles they're just trying to "break things").
+    if w.bit_length() >= 46:
+        # "Only" had < 53 - 46 = 7 bits to spare in IEEE-754 double.
+        # XXX I can't test this - don't have 169 terabytes of RAM to
+        # build a string long enough to trigger this.
+        raise ValueError(f"cannot convert string of len {len(s)} to int")
+    with decimal.localcontext(_unbounded_dec_context) as ctx:
+        D256 = D(256)
+        pow256 = compute_powers(w, D256, BYTELIM)
+        rpow256 = compute_powers(w, 1 / D256, BYTELIM)
+        # We're going to do inexact, chopped arithmetic, multiplying by
+        # an approximation to the reciprocal of 256**i. We chop to get a
+        # lower bound on the true integer quotient. Our approximation is
+        # a lower bound, the multiply is chopped too, and
+        # to_integral_value() is also chopped.
+        ctx.traps[decimal.Inexact] = 0
+        ctx.rounding = decimal.ROUND_DOWN
+        for k, v in pow256.items():
+            # No need to save more precision in the reciprocal than the
+            # power of 256 has, plus some guard digits to absorb most
+            # relevant rounding errors. This is highly signficant:
+            # 1/2**i has the same number of significant decimal digits
+            # as 5**i, generally over twice the number in 2**i,
+            ctx.prec = v.adjusted() + 8
+            # The unary "+" chope the reciprocal back to that precision.
+            pow256[k] = v, +rpow256[k]
+        del rpow256 # exact reciprocals no longer needed
+        ctx.prec = decimal.MAX_PREC
+        inner(D(s), w)
+    return int.from_bytes(result)
 
 def int_from_string(s):
     """Asymptotically fast version of PyLong_FromString(), conversion
@@ -219,7 +334,10 @@ def int_from_string(s):
     # and underscores, and stripped leading whitespace.  The input can still
     # contain underscores and have trailing whitespace.
     s = s.rstrip().replace('_', '')
-    return _str_to_int_inner(s)
+    func = _str_to_int_inner
+    if len(s) >= 3_500_000 and _decimal is not None:
+        func = _dec_str_to_int_inner
+    return func(s)
 
 def str_to_int(s):
     """Asymptotically fast version of decimal string to 'int' conversion."""
@@ -361,3 +479,21 @@ def int_divmod(a, b):
         return ~q, b + ~r
     else:
         return _divmod_pos(a, b)
+
+
+# Notes on _dec_str_to_int_inner:
+#
+# Stefan Pochmann worked up a str->int function that used the decimal
+# module to, in effect, convert from base 10 to base 256. This is
+# "unnatural", in that it requires multiplying and dividing by large
+# powers of 2, which `decimal` isn't naturally suited to. But
+# `decimal`'s `*` and `/` are asymptotically superior to CPython's, so
+# at _some_ point it could be expected to win.
+#
+# Alas, the crossover point was too high to be of much real interest. I
+# (Tim) then worked on ways to replace its division with multiplication
+# by a cached reciprocal approximation instead, fixing up errors
+# afterwards. This reduced the crossover point significantly,
+#
+# I revisited tha code, and found ways to improve and simplify it. The
+# crossover point is at about 3.4 million digits now.

Lib/test/test_int.py

@@ -919,5 +919,15 @@ def test_pylong_roundtrip(self):
             self.assertEqual(n, int(sn))
             bits <<= 1
 
+    @support.requires_resource('cpu')
+    def test_pylong_roundtrip_huge(self):
+        # k blocks of 1234567890
+        k = 1_000_000 # so 10 million digits in all
+        tentoten = 10**10
+        n = 1234567890 * ((tentoten**k - 1) // (tentoten - 1))
+        sn = "1234567890" * k
+        self.assertEqual(n, int(sn))
+        self.assertEqual(sn, str(n))
+
 if __name__ == "__main__":
     unittest.main()

Misc/NEWS.d/next/Core and Builtins/2024-05-09-02-37-25.gh-issue-118750.7aLfT-.rst

@@ -0,0 +1 @@
+If the C version of the ``decimal`` module is available, ``int(str)`` now uses it to supply an asymptotically much faster conversion. However, this only applies if the string contains over about 3.5 million digits.