mirror of https://github.com/python/cpython.git
730 lines
28 KiB
Python
730 lines
28 KiB
Python
"""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 re
|
|
import decimal
|
|
try:
|
|
import _decimal
|
|
except ImportError:
|
|
_decimal = None
|
|
|
|
# A number of functions have this form, where `w` is a desired number of
|
|
# digits in base `base`:
|
|
#
|
|
# def inner(...w...):
|
|
# if w <= LIMIT:
|
|
# return something
|
|
# lo = w >> 1
|
|
# hi = w - lo
|
|
# something involving base**lo, inner(...lo...), j, and inner(...hi...)
|
|
# figure out largest w needed
|
|
# result = inner(w)
|
|
#
|
|
# They all had some on-the-fly scheme to cache `base**lo` results for reuse.
|
|
# Power is costly.
|
|
#
|
|
# This routine aims to compute all amd only the needed powers in advance, as
|
|
# efficiently as reasonably possible. This isn't trivial, and all the
|
|
# on-the-fly methods did needless work in many cases. The driving code above
|
|
# changes to:
|
|
#
|
|
# figure out largest w needed
|
|
# mycache = compute_powers(w, base, LIMIT)
|
|
# result = inner(w)
|
|
#
|
|
# and `mycache[lo]` replaces `base**lo` in the inner function.
|
|
#
|
|
# If an algorithm wants the powers of ceiling(w/2) instead of the floor,
|
|
# pass keyword argument `need_hi=True`.
|
|
#
|
|
# While this does give minor speedups (a few percent at best), the
|
|
# primary intent is to simplify the functions using this, by eliminating
|
|
# the need for them to craft their own ad-hoc caching schemes.
|
|
#
|
|
# See code near end of file for a block of code that can be enabled to
|
|
# run millions of tests.
|
|
def compute_powers(w, base, more_than, *, need_hi=False, show=False):
|
|
seen = set()
|
|
need = set()
|
|
ws = {w}
|
|
while ws:
|
|
w = ws.pop() # any element is fine to use next
|
|
if w in seen or w <= more_than:
|
|
continue
|
|
seen.add(w)
|
|
lo = w >> 1
|
|
hi = w - lo
|
|
# only _need_ one here; the other may, or may not, be needed
|
|
which = hi if need_hi else lo
|
|
need.add(which)
|
|
ws.add(which)
|
|
if lo != hi:
|
|
ws.add(w - which)
|
|
|
|
# `need` is the set of exponents needed. To compute them all
|
|
# efficiently, possibly add other exponents to `extra`. The goal is
|
|
# to ensure that each exponent can be gotten from a smaller one via
|
|
# multiplying by the base, squaring it, or squaring and then
|
|
# multiplying by the base.
|
|
#
|
|
# If need_hi is False, this is already the case (w can always be
|
|
# gotten from w >> 1 via one of the squaring strategies). But we do
|
|
# the work anyway, just in case ;-)
|
|
#
|
|
# Note that speed is irrelevant. These loops are working on little
|
|
# ints (exponents) and go around O(log w) times. The total cost is
|
|
# insignificant compared to just one of the bigint multiplies.
|
|
cands = need.copy()
|
|
extra = set()
|
|
while cands:
|
|
w = max(cands)
|
|
cands.remove(w)
|
|
lo = w >> 1
|
|
if lo > more_than and w-1 not in cands and lo not in cands:
|
|
extra.add(lo)
|
|
cands.add(lo)
|
|
assert need_hi or not extra
|
|
|
|
d = {}
|
|
for n in sorted(need | extra):
|
|
lo = n >> 1
|
|
hi = n - lo
|
|
if n-1 in d:
|
|
if show:
|
|
print("* base", end="")
|
|
result = d[n-1] * base # cheap!
|
|
elif lo in d:
|
|
# Multiplying a bigint by itself is about twice as fast
|
|
# in CPython provided it's the same object.
|
|
if show:
|
|
print("square", end="")
|
|
result = d[lo] * d[lo] # same object
|
|
if hi != lo:
|
|
if show:
|
|
print(" * base", end="")
|
|
assert 2 * lo + 1 == n
|
|
result *= base
|
|
else: # rare
|
|
if show:
|
|
print("pow", end='')
|
|
result = base ** n
|
|
if show:
|
|
print(" at", n, "needed" if n in need else "extra")
|
|
d[n] = result
|
|
|
|
assert need <= d.keys()
|
|
if excess := d.keys() - need:
|
|
assert need_hi
|
|
for n in excess:
|
|
del d[n]
|
|
return d
|
|
|
|
_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
|
|
|
|
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_.
|
|
|
|
from decimal import Decimal as D
|
|
BITLIM = 200
|
|
|
|
# Don't bother caching the "lo" mask in this; the time to compute it is
|
|
# tiny compared to the multiply.
|
|
def inner(n, w):
|
|
if w <= BITLIM:
|
|
return D(n)
|
|
w2 = w >> 1
|
|
hi = n >> w2
|
|
lo = n & ((1 << w2) - 1)
|
|
return inner(lo, w2) + inner(hi, w - w2) * w2pow[w2]
|
|
|
|
with decimal.localcontext(_unbounded_dec_context):
|
|
nbits = n.bit_length()
|
|
w2pow = compute_powers(nbits, D(2), BITLIM)
|
|
if n < 0:
|
|
negate = True
|
|
n = -n
|
|
else:
|
|
negate = False
|
|
result = inner(n, nbits)
|
|
if negate:
|
|
result = -result
|
|
return result
|
|
|
|
def int_to_decimal_string(n):
|
|
"""Asymptotically fast conversion of an 'int' to a decimal string."""
|
|
w = n.bit_length()
|
|
if w > 450_000 and _decimal is not None:
|
|
# It is only usable with the C decimal implementation.
|
|
# _pydecimal.py calls str() on very large integers, which in its
|
|
# turn calls int_to_decimal_string(), causing very deep recursion.
|
|
return str(int_to_decimal(n))
|
|
|
|
# Fallback algorithm for the case when the C decimal module isn't
|
|
# available. This algorithm is asymptotically worse than the algorithm
|
|
# using the decimal module, but better than the quadratic time
|
|
# implementation in longobject.c.
|
|
|
|
DIGLIM = 1000
|
|
def inner(n, w):
|
|
if w <= DIGLIM:
|
|
return str(n)
|
|
w2 = w >> 1
|
|
hi, lo = divmod(n, pow10[w2])
|
|
return inner(hi, w - w2) + inner(lo, w2).zfill(w2)
|
|
|
|
# The estimation of the number of decimal digits.
|
|
# There is no harm in small error. If we guess too large, there may
|
|
# be leading 0's that need to be stripped. If we guess too small, we
|
|
# may need to call str() recursively for the remaining highest digits,
|
|
# which can still potentially be a large integer. This is manifested
|
|
# only if the number has way more than 10**15 digits, that exceeds
|
|
# the 52-bit physical address limit in both Intel64 and AMD64.
|
|
w = int(w * 0.3010299956639812 + 1) # log10(2)
|
|
pow10 = compute_powers(w, 5, DIGLIM)
|
|
for k, v in pow10.items():
|
|
pow10[k] = v << k # 5**k << k == 5**k * 2**k == 10**k
|
|
if n < 0:
|
|
n = -n
|
|
sign = '-'
|
|
else:
|
|
sign = ''
|
|
s = inner(n, w)
|
|
if s[0] == '0' and n:
|
|
# If our guess of w is too large, there may be leading 0's that
|
|
# need to be stripped.
|
|
s = s.lstrip('0')
|
|
return sign + s
|
|
|
|
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
|
|
|
|
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)))
|
|
|
|
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
|
|
|
|
# _spread is for internal testing. It maps a key to the number of times
|
|
# that condition obtained in _dec_str_to_int_inner:
|
|
# key 0 - quotient guess was right
|
|
# key 1 - quotient had to be boosted by 1, one time
|
|
# key 999 - one adjustment wasn't enough, so fell back to divmod
|
|
from collections import defaultdict
|
|
_spread = defaultdict(int)
|
|
del defaultdict
|
|
|
|
def _dec_str_to_int_inner(s, *, GUARD=8):
|
|
# Yes, BYTELIM is "large". Large enough that CPython will usually
|
|
# use the Karatsuba _str_to_int_inner to convert the string. This
|
|
# allowed reducing the cutoff for calling _this_ function from 3.5M
|
|
# to 2M digits. We could almost certainly do even better by
|
|
# fine-tuning this and/or using a larger output base than 256.
|
|
BYTELIM = 100_000
|
|
D = decimal.Decimal
|
|
result = bytearray()
|
|
# See notes at end of file for discussion of GUARD.
|
|
assert GUARD > 0 # if 0, `decimal` can blow up - .prec 0 not allowed
|
|
|
|
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))`.
|
|
# Worse, `int(Decimal) is still quadratic-time for much
|
|
# larger ints. So unless/until all that is repaired, the
|
|
# seemingly redundant `str(Decimal)` is crucial to speed.
|
|
result.extend(int(str(n)).to_bytes(w)) # big-endian default
|
|
return
|
|
w1 = w >> 1
|
|
w2 = w - w1
|
|
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,
|
|
ctx.prec = max(n.adjusted() - p256.adjusted(), 0) + GUARD
|
|
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 on inputs through 5 billion digit strings, the
|
|
# test triggered once in about 200 thousand tries.
|
|
count = 0
|
|
if lo >= p256:
|
|
count = 1
|
|
lo -= p256
|
|
hi += 1
|
|
if lo >= p256:
|
|
# Complete correction via an exact computation. I
|
|
# believe it's not possible to get here provided
|
|
# GUARD >= 3. It's tested by reducing GUARD below
|
|
# that.
|
|
count = 999
|
|
hi2, lo = divmod(lo, p256)
|
|
hi += hi2
|
|
_spread[count] += 1
|
|
# The assert should always succeed, but way too slow to keep
|
|
# enabled.
|
|
#assert hi, lo == divmod(n, pow256[w2][0])
|
|
inner(hi, w1)
|
|
del hi # at top levels, can free a lot of RAM "early"
|
|
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 bound, 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.
|
|
|
|
# The "WASI" test platform can complain about `len(s)` if it's too
|
|
# large to fit in its idea of "an index-sized integer".
|
|
lenS = s.__len__()
|
|
log_ub = lenS * _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've 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 (unless they're just trying to "break things").
|
|
if w.bit_length() >= 46:
|
|
# "Only" had < 53 - 46 = 7 bits to spare in IEEE-754 double.
|
|
raise ValueError(f"cannot convert string of len {lenS} to int")
|
|
with decimal.localcontext(_unbounded_dec_context) as ctx:
|
|
D256 = D(256)
|
|
pow256 = compute_powers(w, D256, BYTELIM, need_hi=True)
|
|
rpow256 = compute_powers(w, 1 / D256, BYTELIM, need_hi=True)
|
|
# 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 multiplication 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 much more precision in the reciprocal than
|
|
# the power of 256 has, plus some guard digits to absorb
|
|
# most relevant rounding errors. This is highly significant:
|
|
# 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() + GUARD + 1
|
|
# The unary "+" chops 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
|
|
of a string of decimal digits into an 'int'."""
|
|
# 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('_', '')
|
|
func = _str_to_int_inner
|
|
if len(s) >= 2_000_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."""
|
|
# 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 b == 0:
|
|
raise ZeroDivisionError('division by zero')
|
|
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)
|
|
|
|
|
|
# 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 the code, and found ways to improve and simplify it. The
|
|
# crossover point is at about 3.4 million digits now.
|
|
#
|
|
# About .adjusted()
|
|
# -----------------
|
|
# Restrict to Decimal values x > 0. We don't use negative numbers in the
|
|
# code, and I don't want to have to keep typing, e.g., "absolute value".
|
|
#
|
|
# For convenience, I'll use `x.a` to mean `x.adjusted()`. x.a doesn't
|
|
# look at the digits of x, but instead returns an integer giving x's
|
|
# order of magnitude. These are equivalent:
|
|
#
|
|
# - x.a is the power-of-10 exponent of x's most significant digit.
|
|
# - x.a = the infinitely precise floor(log10(x))
|
|
# - x can be written in this form, where f is a real with 1 <= f < 10:
|
|
# x = f * 10**x.a
|
|
#
|
|
# Observation; if x is an integer, len(str(x)) = x.a + 1.
|
|
#
|
|
# Lemma 1: (x * y).a = x.a + y.a, or one larger
|
|
#
|
|
# Proof: Write x = f * 10**x.a and y = g * 10**y.a, where f and g are in
|
|
# [1, 10). Then x*y = f*g * 10**(x.a + y.a), where 1 <= f*g < 100. If
|
|
# f*g < 10, (x*y).a is x.a+y.a. Else divide f*g by 10 to bring it back
|
|
# into [1, 10], and add 1 to the exponent to compensate. Then (x*y).a is
|
|
# x.a+y.a+1.
|
|
#
|
|
# Lemma 2: ceiling(log10(x/y)) <= x.a - y.a + 1
|
|
#
|
|
# Proof: Express x and y as in Lemma 1. Then x/y = f/g * 10**(x.a -
|
|
# y.a), where 1/10 < f/g < 10. If 1 <= f/g, (x/y).a is x.a-y.a. Else
|
|
# multiply f/g by 10 to bring it back into [1, 10], and subtract 1 from
|
|
# the exponent to compensate. Then (x/y).a is x.a-y.a-1. So the largest
|
|
# (x/y).a can be is x.a-y.a. Since that's the floor of log10(x/y). the
|
|
# ceiling is at most 1 larger (with equality iff f/g = 1 exactly).
|
|
#
|
|
# GUARD digits
|
|
# ------------
|
|
# We only want the integer part of divisions, so don't need to build
|
|
# the full multiplication tree. But using _just_ the number of
|
|
# digits expected in the integer part ignores too much. What's left
|
|
# out can have a very significant effect on the quotient. So we use
|
|
# GUARD additional digits.
|
|
#
|
|
# The default 8 is more than enough so no more than 1 correction step
|
|
# was ever needed for all inputs tried through 2.5 billion digits. In
|
|
# fact, I believe 3 guard digits are always enough - but the proof is
|
|
# very involved, so better safe than sorry.
|
|
#
|
|
# Short course:
|
|
#
|
|
# If prec is the decimal precision in effect, and we're rounding down,
|
|
# the result of an operation is exactly equal to the infinitely precise
|
|
# result times 1-e for some real e with 0 <= e < 10**(1-prec). In
|
|
#
|
|
# ctx.prec = max(n.adjusted() - p256.adjusted(), 0) + GUARD
|
|
# hi = +n * +recip # unary `+` chops to ctx.prec digits
|
|
#
|
|
# we have 3 visible chopped operations, but there's also a 4th:
|
|
# precomputing a truncated `recip` as part of setup.
|
|
#
|
|
# So the computed product is exactly equal to the true product times
|
|
# (1-e1)*(1-e2)*(1-e3)*(1-e4); since the e's are all very small, an
|
|
# excellent approximation to the second factor is 1-(e1+e2+e3+e4) (the
|
|
# 2nd and higher order terms in the expanded product are too tiny to
|
|
# matter). If they're all as large as possible, that's
|
|
#
|
|
# 1 - 4*10**(1-prec). This, BTW, is all bog-standard FP error analysis.
|
|
#
|
|
# That implies the computed product is within 1 of the true product
|
|
# provided prec >= log10(true_product) + 1.602.
|
|
#
|
|
# Here are telegraphic details, rephrasing the initial condition in
|
|
# equivalent ways, step by step:
|
|
#
|
|
# prod - prod * (1 - 4*10**(1-prec)) <= 1
|
|
# prod - prod + prod * 4*10**(1-prec)) <= 1
|
|
# prod * 4*10**(1-prec)) <= 1
|
|
# 10**(log10(prod)) * 4*10**(1-prec)) <= 1
|
|
# 4*10**(1-prec+log10(prod))) <= 1
|
|
# 10**(1-prec+log10(prod))) <= 1/4
|
|
# 1-prec+log10(prod) <= log10(1/4) = -0.602
|
|
# -prec <= -1.602 - log10(prod)
|
|
# prec >= log10(prod) + 1.602
|
|
#
|
|
# The true product is the same as the true ratio n/p256. By Lemma 2
|
|
# above, n.a - p256.a + 1 is an upper bound on the ceiling of
|
|
# log10(prod). Then 2 is the ceiling of 1.602. so n.a - p256.a + 3 is an
|
|
# upper bound on the right hand side of the inequality. Any prec >= that
|
|
# will work.
|
|
#
|
|
# But since this is just a sketch of a proof ;-), the code uses the
|
|
# empirically tested 8 instead of 3. 5 digits more or less makes no
|
|
# practical difference to speed - these ints are huge. And while
|
|
# increasing GUARD above 3 may not be necessary, every increase cuts the
|
|
# percentage of cases that need a correction at all.
|
|
#
|
|
# On Computing Reciprocals
|
|
# ------------------------
|
|
# In general, the exact reciprocals we compute have over twice as many
|
|
# significant digits as needed. 1/256**i has the same number of
|
|
# significant decimal digits as 5**i. It's a significant waste of RAM
|
|
# to store all those unneeded digits.
|
|
#
|
|
# So we cut exact reciprocals back to the least precision that can
|
|
# be needed so that the error analysis above is valid,
|
|
#
|
|
# [Note: turns out it's very significantly faster to do it this way than
|
|
# to compute 1 / 256**i directly to the desired precision, because the
|
|
# power method doesn't require division. It's also faster than computing
|
|
# (1/256)**i directly to the desired precision - no material division
|
|
# there, but `compute_powers()` is much smarter about _how_ to compute
|
|
# all the powers needed than repeated applications of `**` - that
|
|
# function invokes `**` for at most the few smallest powers needed.]
|
|
#
|
|
# The hard part is that chopping back to a shorter width occurs
|
|
# _outside_ of `inner`. We can't know then what `prec` `inner()` will
|
|
# need. We have to pick, for each value of `w2`, the largest possible
|
|
# value `prec` can become when `inner()` is working on `w2`.
|
|
#
|
|
# This is the `prec` inner() uses:
|
|
# max(n.a - p256.a, 0) + GUARD
|
|
# and what setup uses (renaming its `v` to `p256` - same thing):
|
|
# p256.a + GUARD + 1
|
|
#
|
|
# We need that the second is always at least as large as the first,
|
|
# which is the same as requiring
|
|
#
|
|
# n.a - 2 * p256.a <= 1
|
|
#
|
|
# What's the largest n can be? n < 255**w = 256**(w2 + (w - w2)). The
|
|
# worst case in this context is when w ix even. and then w = 2*w2, so
|
|
# n < 256**(2*w2) = (256**w2)**2 = p256**2. By Lemma 1, then, n.a
|
|
# is at most p256.a + p256.a + 1.
|
|
#
|
|
# So the most n.a - 2 * p256.a can be is
|
|
# p256.a + p256.a + 1 - 2 * p256.a = 1. QED
|
|
#
|
|
# Note: an earlier version of the code split on floor(e/2) instead of on
|
|
# the ceiling. The worst case then is odd `w`, and a more involved proof
|
|
# was needed to show that adding 4 (instead of 1) may be necessary.
|
|
# Basically because, in that case, n may be up to 256 times larger than
|
|
# p256**2. Curiously enough, by splitting on the ceiling instead,
|
|
# nothing in any proof here actually depends on the output base (256).
|
|
|
|
# Enable for brute-force testing of compute_powers(). This takes about a
|
|
# minute, because it tries millions of cases.
|
|
if 0:
|
|
def consumer(w, limit, need_hi):
|
|
seen = set()
|
|
need = set()
|
|
def inner(w):
|
|
if w <= limit:
|
|
return
|
|
if w in seen:
|
|
return
|
|
seen.add(w)
|
|
lo = w >> 1
|
|
hi = w - lo
|
|
need.add(hi if need_hi else lo)
|
|
inner(lo)
|
|
inner(hi)
|
|
inner(w)
|
|
exp = compute_powers(w, 1, limit, need_hi=need_hi)
|
|
assert exp.keys() == need
|
|
|
|
from itertools import chain
|
|
for need_hi in (False, True):
|
|
for limit in (0, 1, 10, 100, 1_000, 10_000, 100_000):
|
|
for w in chain(range(1, 100_000),
|
|
(10**i for i in range(5, 30))):
|
|
consumer(w, limit, need_hi)
|