highest power of of prime number $p$ among $k$ consecutive positive integers

Question

What is a good way to determine the highest power of of prime number $p$ among $k$ consecutive positive integers $n$ to $n+k$?

For example:

• consecutive integers $\{12,13,14,15,16,17,18,19\}$

• highest power of $2$: $4$ at $16$

• highest power of $3$: $2$ at $18$

Answer 1

You can try with a simple recursive algorithm: once you know that $p$ divides one of the numbers you immediately know which numbers are divided by $p$. If you discard the other numbers you can then iterate this process until you are left with only one number, so keep dividing it by $p$ until it is not divisible by it anymore. Note that this doesn't assume $p$ to be prime.

Here's a quick and dirty Python script implementing this algorithm. It can likely be improved, but so far it seems to work well enough.

"""
Usage: highest_power.py <d> <n> <k>

Returns the time taken by 100 iterations of highest_power(<d>, <n>, <k>).
"""

def divides(d,n):
    return n % d == 0

def highest_power(d,n,k):
    """
    Returns a tuple '(p,i)' where 'p' is the highest power of 'd' between
    'n' and the successive 'k' integers, while 'i' is the smallest of
    those integers divided by 'd ** p'.
    """

    integers = range(n,n+k)
    power = 0

    while True:
        minimum = integers[0]
        integers = [x//d for x in integers if divides(d,x)]

        if integers == []:
            return (power, minimum * (d ** power))
        else:
            power += 1

if __name__ == "__main__":
    from sys import argv
    import timeit

    script, d, n, k = argv

    print(timeit.timeit("highest_power({0}, {1}, {2})".format(d, n, k),
                        setup="from __main__ import highest_power",
                        number=100))

And here's a simple benchmark:

$python highest_power_1.py 2 1 $(python -c 'print(10)')
0.00114600800043263

$python highest_power_1.py 2 1 $(python -c 'print(10**3)')
0.05992611500005296

$python highest_power_1.py 2 1 $(python -c 'print(10**6)')
35.90007105699988

Remark: I think that the main bottleneck of the above code is the use of lists: for $k = 10^9$ this consumes way too much memory to be usable on my system. I don't know if there is a clean way to avoid them, though.

---

Update: A more efficient way is to recursively enumerate, starting at $n$, the smallest numbers divisible by the successive powers of $p$, stopping as soon as we reach $n+k$. Mathematically this is fairly easy: if $d$ and $m$ are two given positive integers such that $d \nmid m$, then $$ m + d - (m \mod d) $$ is the first number after $m$ divisible by $d$.

Here's a Python script implementing this algorithm -- note that I (ab)used generators to simplify the book-keeping and to split the code in small, easily testable pieces.

"""
Usage: highest_power.py <d> <n> <k>

Returns the time taken by 100 iterations of highest_power(<d>, <n>, <k>).
"""

from more_itertools import peekable, iterate


def take_one_while(predicate, iterable):
    """
    Return the last element from the head of 'iterable' for which
    'predicate' is true. If 'predicate' is false for the first element
    of 'iterable' return 'None'. Equivalent to:

    from itertools import takewhile

    def take_one_while(predicate, iterable):
        # take_one_while(lambda x: x<5, [1,2,4,6,4,1]) --> 4
        # take_one_while(lambda x: x<1, [2,3,0]) --> None

        try:
            return list(takewhile(predicate, iterable))[-1]
        except IndexError:
            return None
    """
    it = peekable(iterable)
    result = None

    while predicate(it.peek()):
        result = next(it)

    return result


def max_factor(d, n, start=1):
    """Return the maximum power of 'd' that divides 'n'.

    If 'start' is given, return 'start * (d**s)' where 's' is such that
    'start * (d ** s)' divides 'n' but 'start * (d ** (s + 1))' does not.
    Return 'None' if 'start' doesn't divide 'n'.
    """
    powers = iterate(lambda x: x*d, start)

    def divides_n(p):
        return n % p == 0

    return take_one_while(divides_n, powers)


def next_highest_powers(d, start=1):
    """
    Generate the smallest numbers divisible by successive powers of
    'd', beginning at 'start', and the powers of 'd' dividing them.

    'start' defaults to '1', effectively generating the successive
    powers of 'd'.
    """
    power = max_factor(d, start)
    if start % d == 0:
        yield (start, power)

    while True:
        next_power = power * d
        remainder = start % next_power
        start += next_power - remainder
        power = max_factor(d, start)
        yield (start, power)


def highest_power(d, n, k):
    """
    Return '(h, p)', where 'h' is the smallest number between 'n' and
    'n+k-1' divisible by the highest power of 'd'.
    If no number between 'n' and 'n+k-1' is divisible by 'd', return 'None'.
    """
    it = next_highest_powers(d, n)
    return take_one_while(lambda x: x[0] < n+k, it)


if __name__ == "__main__":
    from sys import argv
    import timeit

    script, d, n, k = argv

    print(timeit.timeit("highest_power({0}, {1}, {2})".format(d, n, k),
                        setup="from __main__ import highest_power",
                        number=100))

I'm not sure about how to compute its complexity, but this piece of code is definitely more efficient than the previous one (both time and space wise).

$python highest_power.py 2 1 $(python -c 'print(10)')      
0.007290727000508923

$python highest_power.py 2 1 $(python -c 'print(10**3)')
0.027383217000533477

$python highest_power.py 2 1 $(python -c 'print(10**9)')
0.15067375299986452

$python highest_power.py 2 1 $(python -c 'print(10**12)')
0.2264490380002826

$python highest_power.py 2 1 $(python -c 'print(10**15)')
0.3507304409995413

$python highest_power.py 2 1 $(python -c 'print(10**100)')
12.48393625099925