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