gh-90716: add _pylong.py module (#96673)

Add Python implementations of certain longobject.c functions. These use
asymptotically faster algorithms that can be used for operations on
integers with many digits. In those cases, the performance overhead of
the Python implementation is not significant since the asymptotic
behavior is what dominates runtime. Functions provided by this module
should be considered private and not part of any public API.

Co-author: Tim Peters <tim.peters@gmail.com>
Co-author: Mark Dickinson <dickinsm@gmail.com>
Co-author: Bjorn Martinsson

1 files changed, 295 insertions, 0 deletions

diff --git a/Lib/_pylong.py b/Lib/_pylong.py
new file mode 100644
index 0000000..e0d4e90
--- /dev/null
+++ b/Lib/_pylong.py
@@ -0,0 +1,295 @@
+"""Python implementations of some algorithms for use by longobject.c.
+The goal is to provide asymptotically faster algorithms that can be
+used for operations on integers with many digits.  In those cases, the
+performance overhead of the Python implementation is not significant
+since the asymptotic behavior is what dominates runtime. Functions
+provided by this module should be considered private and not part of any
+public API.
+
+Note: for ease of maintainability, please prefer clear code and avoid
+"micro-optimizations".  This module will only be imported and used for
+integers with a huge number of digits.  Saving a few microseconds with
+tricky or non-obvious code is not worth it.  For people looking for
+maximum performance, they should use something like gmpy2."""
+
+import sys
+import re
+import decimal
+
+_DEBUG = False
+
+
+def int_to_decimal(n):
+    """Asymptotically fast conversion of an 'int' to Decimal."""
+
+    # Function due to Tim Peters.  See GH issue #90716 for details.
+    # https://github.com/python/cpython/issues/90716
+    #
+    # The implementation in longobject.c of base conversion algorithms
+    # between power-of-2 and non-power-of-2 bases are quadratic time.
+    # This function implements a divide-and-conquer algorithm that is
+    # faster for large numbers.  Builds an equal decimal.Decimal in a
+    # "clever" recursive way.  If we want a string representation, we
+    # apply str to _that_.
+
+    if _DEBUG:
+        print('int_to_decimal', n.bit_length(), file=sys.stderr)
+
+    D = decimal.Decimal
+    D2 = D(2)
+
+    BITLIM = 128
+
+    mem = {}
+
+    def w2pow(w):
+        """Return D(2)**w and store the result. Also possibly save some
+        intermediate results. In context, these are likely to be reused
+        across various levels of the conversion to Decimal."""
+        if (result := mem.get(w)) is None:
+            if w <= BITLIM:
+                result = D2**w
+            elif w - 1 in mem:
+                result = (t := mem[w - 1]) + t
+            else:
+                w2 = w >> 1
+                # If w happens to be odd, w-w2 is one larger then w2
+                # now. Recurse on the smaller first (w2), so that it's
+                # in the cache and the larger (w-w2) can be handled by
+                # the cheaper `w-1 in mem` branch instead.
+                result = w2pow(w2) * w2pow(w - w2)
+            mem[w] = result
+        return result
+
+    def inner(n, w):
+        if w <= BITLIM:
+            return D(n)
+        w2 = w >> 1
+        hi = n >> w2
+        lo = n - (hi << w2)
+        return inner(lo, w2) + inner(hi, w - w2) * w2pow(w2)
+
+    with decimal.localcontext() as ctx:
+        ctx.prec = decimal.MAX_PREC
+        ctx.Emax = decimal.MAX_EMAX
+        ctx.Emin = decimal.MIN_EMIN
+        ctx.traps[decimal.Inexact] = 1
+
+        if n < 0:
+            negate = True
+            n = -n
+        else:
+            negate = False
+        result = inner(n, n.bit_length())
+        if negate:
+            result = -result
+    return result
+
+
+def int_to_decimal_string(n):
+    """Asymptotically fast conversion of an 'int' to a decimal string."""
+    return str(int_to_decimal(n))
+
+
+def _str_to_int_inner(s):
+    """Asymptotically fast conversion of a 'str' to an 'int'."""
+
+    # Function due to Bjorn Martinsson.  See GH issue #90716 for details.
+    # https://github.com/python/cpython/issues/90716
+    #
+    # The implementation in longobject.c of base conversion algorithms
+    # between power-of-2 and non-power-of-2 bases are quadratic time.
+    # This function implements a divide-and-conquer algorithm making use
+    # of Python's built in big int multiplication. Since Python uses the
+    # Karatsuba algorithm for multiplication, the time complexity
+    # of this function is O(len(s)**1.58).
+
+    DIGLIM = 2048
+
+    mem = {}
+
+    def w5pow(w):
+        """Return 5**w and store the result.
+        Also possibly save some intermediate results. In context, these
+        are likely to be reused across various levels of the conversion
+        to 'int'.
+        """
+        if (result := mem.get(w)) is None:
+            if w <= DIGLIM:
+                result = 5**w
+            elif w - 1 in mem:
+                result = mem[w - 1] * 5
+            else:
+                w2 = w >> 1
+                # If w happens to be odd, w-w2 is one larger then w2
+                # now. Recurse on the smaller first (w2), so that it's
+                # in the cache and the larger (w-w2) can be handled by
+                # the cheaper `w-1 in mem` branch instead.
+                result = w5pow(w2) * w5pow(w - w2)
+            mem[w] = result
+        return result
+
+    def inner(a, b):
+        if b - a <= DIGLIM:
+            return int(s[a:b])
+        mid = (a + b + 1) >> 1
+        return inner(mid, b) + ((inner(a, mid) * w5pow(b - mid)) << (b - mid))
+
+    return inner(0, len(s))
+
+
+def int_from_string(s):
+    """Asymptotically fast version of PyLong_FromString(), conversion
+    of a string of decimal digits into an 'int'."""
+    if _DEBUG:
+        print('int_from_string', len(s), file=sys.stderr)
+    # PyLong_FromString() has already removed leading +/-, checked for invalid
+    # use of underscore characters, checked that string consists of only digits
+    # 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)
+
+
+def str_to_int(s):
+    """Asymptotically fast version of decimal string to 'int' conversion."""
+    # FIXME: this doesn't support the full syntax that int() supports.
+    m = re.match(r'\s*([+-]?)([0-9_]+)\s*', s)
+    if not m:
+        raise ValueError('invalid literal for int() with base 10')
+    v = int_from_string(m.group(2))
+    if m.group(1) == '-':
+        v = -v
+    return v
+
+
+# Fast integer division, based on code from Mark Dickinson, fast_div.py
+# GH-47701. Additional refinements and optimizations by Bjorn Martinsson.  The
+# algorithm is due to Burnikel and Ziegler, in their paper "Fast Recursive
+# Division".
+
+_DIV_LIMIT = 4000
+
+
+def _div2n1n(a, b, n):
+    """Divide a 2n-bit nonnegative integer a by an n-bit positive integer
+    b, using a recursive divide-and-conquer algorithm.
+
+    Inputs:
+      n is a positive integer
+      b is a positive integer with exactly n bits
+      a is a nonnegative integer such that a < 2**n * b
+
+    Output:
+      (q, r) such that a = b*q+r and 0 <= r < b.
+
+    """
+    if a.bit_length() - n <= _DIV_LIMIT:
+        return divmod(a, b)
+    pad = n & 1
+    if pad:
+        a <<= 1
+        b <<= 1
+        n += 1
+    half_n = n >> 1
+    mask = (1 << half_n) - 1
+    b1, b2 = b >> half_n, b & mask
+    q1, r = _div3n2n(a >> n, (a >> half_n) & mask, b, b1, b2, half_n)
+    q2, r = _div3n2n(r, a & mask, b, b1, b2, half_n)
+    if pad:
+        r >>= 1
+    return q1 << half_n | q2, r
+
+
+def _div3n2n(a12, a3, b, b1, b2, n):
+    """Helper function for _div2n1n; not intended to be called directly."""
+    if a12 >> n == b1:
+        q, r = (1 << n) - 1, a12 - (b1 << n) + b1
+    else:
+        q, r = _div2n1n(a12, b1, n)
+    r = (r << n | a3) - q * b2
+    while r < 0:
+        q -= 1
+        r += b
+    return q, r
+
+
+def _int2digits(a, n):
+    """Decompose non-negative int a into base 2**n
+
+    Input:
+      a is a non-negative integer
+
+    Output:
+      List of the digits of a in base 2**n in little-endian order,
+      meaning the most significant digit is last. The most
+      significant digit is guaranteed to be non-zero.
+      If a is 0 then the output is an empty list.
+
+    """
+    a_digits = [0] * ((a.bit_length() + n - 1) // n)
+
+    def inner(x, L, R):
+        if L + 1 == R:
+            a_digits[L] = x
+            return
+        mid = (L + R) >> 1
+        shift = (mid - L) * n
+        upper = x >> shift
+        lower = x ^ (upper << shift)
+        inner(lower, L, mid)
+        inner(upper, mid, R)
+
+    if a:
+        inner(a, 0, len(a_digits))
+    return a_digits
+
+
+def _digits2int(digits, n):
+    """Combine base-2**n digits into an int. This function is the
+    inverse of `_int2digits`. For more details, see _int2digits.
+    """
+
+    def inner(L, R):
+        if L + 1 == R:
+            return digits[L]
+        mid = (L + R) >> 1
+        shift = (mid - L) * n
+        return (inner(mid, R) << shift) + inner(L, mid)
+
+    return inner(0, len(digits)) if digits else 0
+
+
+def _divmod_pos(a, b):
+    """Divide a non-negative integer a by a positive integer b, giving
+    quotient and remainder."""
+    # Use grade-school algorithm in base 2**n, n = nbits(b)
+    n = b.bit_length()
+    a_digits = _int2digits(a, n)
+
+    r = 0
+    q_digits = []
+    for a_digit in reversed(a_digits):
+        q_digit, r = _div2n1n((r << n) + a_digit, b, n)
+        q_digits.append(q_digit)
+    q_digits.reverse()
+    q = _digits2int(q_digits, n)
+    return q, r
+
+
+def int_divmod(a, b):
+    """Asymptotically fast replacement for divmod, for 'int'.
+    Its time complexity is O(n**1.58), where n = #bits(a) + #bits(b).
+    """
+    if _DEBUG:
+        print('int_divmod', a.bit_length(), b.bit_length(), file=sys.stderr)
+    if b == 0:
+        raise ZeroDivisionError
+    elif b < 0:
+        q, r = int_divmod(-a, -b)
+        return q, -r
+    elif a < 0:
+        q, r = int_divmod(~a, b)
+        return ~q, b + ~r
+    else:
+        return _divmod_pos(a, b)