4
\$\begingroup\$

Background

(Feel free to skip; the help center just recommended adding background information :)

A few days ago, I was looking at some code which happened to use a flat array to store a complete binary tree. That got me thinking about calculating some properties of such a tree. The number of nodes in a perfect binary tree of depth n would be 2 ** n - 1, the index of the parent node, given a child node i would be (i - 1) >> 1 and the level of node i in the tree would be (i + 1).bit_length().

What if we wanted to extend the data structure to k-ary trees? The number of nodes in a perfect k-ary tree of depth n would be (k ** n - 1) // (k - 1), the index of the parent node, given a child node i would be (i - 1) // k and the level of node i in the tree would be base_k_digit_length((i + 1) * (k - 1)) (which is sort of the inverse of the number of nodes calculation).

So, what does base_k_digit_length(x) really correspond to? Well, conceptually, it's ceil(log(x + 1, k)).

Integer logarithms

(Essentially, more background information; Again, feel free to skip)

There are 4 different variants of integer logarithm that I consider interesting, which differ on how the real-valued logarithm is rounded to an integer and correspond to the inequality operators >, ≥, <, and ≤ (based on how the integer return value corresponds to the real return value):

  • gt_log(x, k) = ceil(log(x + 1, k))
  • ge_log(x, k) = ceil(log(x, k))
  • lt_log(x, k) = floor(log(x - 1, k))
  • le_log(x, k) = floor(log(x, k))

Unfortunately, because of the inexactness of floating point, those definitions can't be used directly. For example, at least on my machine, log(125, 5) returns 3.0000000000000004 instead of 3, which means such gt_log(124, 5) would return 4 instead of the correct 3.

Now, said functions could be implemented using very simple loops that would have O(log x) running time, but I wanted to opt for an O(1) implementation based on lookup tables indexed on int.bit_length(). (And I did check - int.bit_length() is O(1), since Python knows the length of the buffer used to store the data for an integer, so it only has to find the most significant 1-bit in the most significant element)

Busting out lookup tables might sound silly, since even in Python, the loop-based approach will only take a few microseconds with any realistically sized argument, so while the lookup table based approach is 2-10x faster, it's a few microseconds we're talking about. But I see it this way: If I used loops for this, that could make some O(n) algorithm in the future into a O(n log n) algorithm, even if the difference wouldn't really be much time-wise (ie. other operations would probably usually shadow the cost of the logarithm).

The code

_lut = {}
def _init_next_power(limits, power, value):
        """Insert LUT entries for all bit lengths between the given power and the previous one

        The given power must be exactly one larger than the previous one!
        """
        bitlen = value.bit_length()
        # All numbers (with previous_bitlen < number.bit_length() <= bitlen) less than the given value have an integer logarithm equal to the given power
        while len(limits) <= bitlen:
                limits.append((power, value))
        return bitlen

def init_upto_power(max_power, base):
        """Initialize or extend the range supported by the integer logarithm functions so that all integers up to base ** max_power are supported for the given base"""
        limits = _lut.setdefault(base, [(0, 1), (0, 1)])
        min_power, value = limits[-1]
        for power in range(min_power + 1, max_power + 1):
                value *= base
                _init_next_power(limits, power, value)
        return power, value

def init_upto_bitlen(max_bitlen, base):
        """Initialize or extend the range supported by the integer logarithm functions so that all integers up to the given bit length are supported for the given base"""
        limits = _lut.setdefault(base, [(0, 1), (0, 1)])
        power, value = limits[-1]
        bitlen = len(limits) - 1
        while bitlen < max_bitlen:
                value *= base
                power += 1
                bitlen = _init_next_power(limits, power, value)
        return power, value

def gt_log(value, base):
        """Return the minimum power of base greater than value

        That is, the returned power satisfies: base ** (power - 1) <= value < base ** power
        This is also the exact integer equivalent of: ceil(log(value + 1, base))
        """
        assert value > 0, 'Logarithm is only defined for numbers greater than zero (the power approaches negative infinity as the value approaches zero)'
        power, limit = _lut[base][value.bit_length()]
        return power + (value >= limit)

def ge_log(value, base):
        """Return the minimum power of base greater than or equal to value

        That is, the returned power satisfies: base ** (power - 1) < value <= base ** power
        This is also the exact integer equivalent of: ceil(log(value, base))
        """
        assert value > 0, 'Logarithm is only defined for numbers greater than zero (the power approaches negative infinity as the value approaches zero)'
        power, limit = _lut[base][value.bit_length()]
        return power + (value > limit)

def lt_log(value, base):
        """Return the maximum power of base less than value

        That is, the returned power satisfies: base ** power < value <= base ** (power + 1)
        This is also the exact integer equivalent of: floor(log(value - 1, base))
        """
        assert value > 0, 'Logarithm is only defined for numbers greater than zero (the power approaches negative infinity as the value approaches zero)'
        power, limit = _lut[base][value.bit_length()]
        return power - (value <= limit)

def le_log(value, base):
        """Return the maximum power of base less than or equal to value

        That is, the returned power satisfies: base ** power <= value < base ** (power + 1)
        This is also the exact integer equivalent of: floor(log(value, base))
        """
        assert value > 0, 'Logarithm is only defined for numbers greater than zero (the power approaches negative infinity as the value approaches zero)'
        power, limit = _lut[base][value.bit_length()]
        return power - (value < limit)

if __name__ == '__main__':
        import math

        # When we say ceil(log(value + 1, base)), we really mean ceil(log(value + epsilon, base)), where 0 < epsilon <= 1
        # Since log(0) is undefined, we use epsilon < 1, so we don't get an error in the floor(log(value - epsilon, base)) case
        epsilon = 0.5

        # Since some of the floating point values returned by log(base ** power, base) are not exactly equal to power, we round
        # the return value of log() to the given precision - just enough to get the correct return value with exact powers
        precision = 14

        max_value = 10**6
        for base in range(2, 11):
                print(base)
                init_upto_bitlen(max_value.bit_length(), base)
                for value in range(1, max_value + 1):
                        power = gt_log(value, base)
                        assert base ** (power - 1) <= value < base ** power
                        assert power == math.ceil(round(math.log(value + epsilon, base), precision))

                        power = ge_log(value, base)
                        assert base ** (power - 1) < value <= base ** power
                        assert power == math.ceil(round(math.log(value, base), precision))

                        power = lt_log(value, base)
                        assert base ** power < value <= base ** (power + 1)
                        assert power == math.floor(round(math.log(value - epsilon, base), precision))

                        power = le_log(value, base)
                        assert base ** power <= value < base ** (power + 1)
                        assert power == math.floor(round(math.log(value, base), precision))

I am mostly concerned about the naming of the initialization functions; I feel like init_upto_* sounds like the functions can only be called once, for initialization. Then again, if the names were something like extend_range_to_*, I wonder if it'd be too non-obvious that one of them has to be called for initialization before using the other functions. (I also considered using "by" or "with" in place of "upto" or "to", in the same way you might have group_by_length)

I also thought about some kind of ceil_log() -style naming for the actual lookup functions, but couldn't figure out a reasonable way to differentiate between the > and ≥ cases. The working name for the module is intlog - for "integer logarithm" - but I'm not perfectly convinced about that either.

I fully expect to get comments along the lines of "The xx_log() functions should just dynamically extend the lookup table!", and that is perfectly valid critique, but since the functions are so small, I feel like the percentual slowdown caused by such logic would be a bit regrettable. Still, I might end up making that the default, and adding a fast_ prefix to the current functions. In any case, even if the functions did do the LUT growing automatically, the init functions would still need a good name. (Although the names would be a bit less critical, since they wouldn't be a part of the public API)

While I'm mostly concerned about the naming - since I can see clear problems with it - I'm of course open to any and all possible feedback!

\$\endgroup\$
2
\$\begingroup\$

I think this code is pretty good, I like the docstrings!

Some things strike me as a bit awkward though,

  1. Stay DRY

    These log functions are really similar, the only thing that differentiates is the operators

    You could make the operator an argument of the function, then you'd only have 1 function

    You'd then need a lookup table for the operators so you know when to subtract or add.

  2. The are better testing modules

    Assertions are fine, but are not the best testsuites, there are some good modules for testing in python (I'm a big fan of doctest) but unittest is also pretty good

So I would rewrite the log functions to this:

import operator
import doctest

MATH_OPER = {
                operator.gt: operator.add, 
                operator.ge: operator.add,
                operator.lt: operator.sub,
                operator.le: operator.sub
            }

def log_oper(op, value, base):
    """
    The logarithmic function with an operator

    Test all operators
    >>> log_oper(operator.gt, 3, 2)
    2
    >>> log_oper(operator.ge, 3, 2)
    2
    >>> log_oper(operator.lt, 3, 2)
    1
    >>> log_oper(operator.le, 3, 2)
    1

    Test assertion Exception
    >>> log_oper(operator.le, -1, 2)
    Traceback (most recent call last):
     ...
    AssertionError: Logarithm is only defined for numbers greater than zero (the power approaches negative infinity as the value approaches zero)
    """
    add_or_sub = MATH_OPER[op]
    assert value > 0, 'Logarithm is only defined for numbers greater than zero (the power approaches negative infinity as the value approaches zero)'
    power, limit = _lut[base][value.bit_length()]
    return add_or_sub(power, op(value, limit))

if __name__ == '__main__':
    max_value = 10**6
    base = 2
    init_upto_bitlen(max_value.bit_length(), base)
    doctest.testmod()
\$\endgroup\$
  • \$\begingroup\$ Ah, for some reason it didn't even cross my mind to merge the functions! Although you'd also need to add the comparison operators to the MATH_OPER lookup table, since the operator defining what the log function does is different from the one used in the code. (Ie. for gt_log, you want to return a power greater than the one returned by normal log, which means you have to add one to the power if the given value is greater than or equal to the limit value - and likewise for the other functions; the operator always differs on equality) \$\endgroup\$ – Aleksi Torhamo Sep 22 '18 at 7:41
  • \$\begingroup\$ @AleksiTorhamo Maybe I'm not getting this correctly, but that is what the code does. You can give the operator as a parameter, and op(value, limit) it will be either 1 or 0 depending if it is greater then or not \$\endgroup\$ – Ludisposed Sep 22 '18 at 12:03
  • \$\begingroup\$ As it is, if I had gt_log(value, base), I'd have to replace it with log_oper(operator.ge, value, base) (so gt gets replaced by ge), so the op argument is different from what the function actually does, requiring the user to know and remember the internal details of the implementation to select the correct op argument. So in gt_log(value, base) the gt stands for "Return an integer greater than log(value, base)", whereas in log_oper(operation.ge, value, base), there is no clear relation to the caller between operation.ge and the actual return value of the function. \$\endgroup\$ – Aleksi Torhamo Sep 22 '18 at 12:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.