Project Euler 30 -- Digit fifth powers

Question

Here's my solution for Project Euler Problem 30. Any optimizations in time, space, or style would be greatly appreciated.

from timeit import default_timer as timer

def digit_powers(exponent):
    def power(k):
        return int(k) ** exponent
    if exponent <= 1:
        return "The exponent must be at least 2."
    else:
        total_sum = 0
        upper_bound = (exponent + 1) * (9 ** exponent)
        for number in range(10, upper_bound + 1):
            digits = [x for x in str(number)]
            if number == sum(map(power, digits)):
                total_sum += number
        return total_sum

start = timer()
ans = digit_powers(5)
elapsed_time = (timer() - start) * 1000 # s --> ms

print "Found %d in %r ms." % (ans, elapsed_time)

Answer 1

In looking at your code, it seems well written and solves the problem in a straight forward manner. The one thing I would do is to remove the else after your argument check and remove a level of indent. That just helps to keep things flat... Flat is better than nested.

Now onto optimization and potentially making it more complex than it needs to be in order to make it faster. In these problems it seems like a reasonable thing to do even though....Simple is better than complex.

Here are a couple things I looked at as potential issues:

1. When doing the sum(map(power, digits)), you are calculating the power of each digit many times. If this was a much more expensive calculation then only calculating 1 time is even more important.
2. In this problem there is a conversion from int->str->int. If you can breakdown the digits without converting them it saves some time. It seems like a good thing to know how to do just in case your ever trying a different language where skipping those transitions may be even faster.

To address these, I did two things:

1. Pre-calculate 1-9 powers and store the results in a dictionary so they don't get recalculated (Memoization). It lets me look them up very quickly without calculating each time.
2. Separate the digits using math.log10 as to avoid the int->str->int conversion.

The first item is fairly straight forward, calculating for each digit:

powers = {}
for a in range(10):
    powers[a] = a ** exponent

The second part is a little trickier, what it does is uses math.log10 which gives a count of the digits, then it does a mod 10 to get the right most digit, looks up that power to add to the total and then finally does integer division to divide by 10 effectively moving the decimal place 1 position to the left (dropping the right most digit we just dealt with). This is probably a little different between Python 2/3 due to changes in int/float division. It is basically doing what you are using [x for x in str(number)] to do but without converting to string and adding the power of 5 dict lookup as it goes along:

savei = i
for _ in range(int(math.log10(i)) + 1):
    digit = i % 10
    total += powers[digit]
    i //= 10

Applying those changes to the code make it about twice as fast. Here is what I end up with for the function:

import math

def digit_powers(exponent):
    if exponent <= 1:
        return "The exponent must be at least 2."

    powers = {}
    answer = 0
    # Get the powers
    for a in range(10):
        powers[a] = a ** exponent

    limit = (exponent + 1) * (9 ** exponent)

    # Search for them
    for i in range(10, limit):
        savei = i
        total = 0

        for _ in range(int(math.log10(i)) + 1):
            digit = i % 10
            total += powers[digit]
            i //= 10

        if (total == savei):
            answer += total

    return answer

Answer 2

You can make this much faster by precalculating the powers of each digit. Also note that the order of processing for digits doesn't matter so this code starts with the ones place. Here's what I came up with, and on my machine it ran 4x faster than the original code:

def digit_powers(exponent):      
    if exponent <= 1:
        return "The exponent must be at least 2."
    powdigits = [i**exponent for i in range(10)]   
    total_sum = 0
    upper_bound = (exponent + 1) * powdigits[9]
    for number in range(10, upper_bound + 1):
        partialsum = tempnum = number
        while tempnum:
            partialsum -= powdigits[tempnum%10]
            tempnum /= 10
        if not partialsum:
            total_sum += number
    return total_sum

Answer 3

I doubt you will see my answer. Hopefully it can help someone else in the future. First I want to compare the answers above speedwise. As a tease I will include my own attempt

Comparison

-------------------------------
|  Name           |  Time     |
|  Joshua (OP)    |  2000 ms  |
|  clutton        |   500 ms  |
|  Edward         |   300 ms  |
|  Nebuchadnezzar |    20 ms  |
-------------------------------

So the answers above give a rough 4x-5x speed improvement. So how can we increase this to a 100x speed improvement? Well naively you check every number up to the limit \$ 5*9^5 \$. This means your code does 443839 checks, and for every check it has to compute the digit sum. In order to get a significant speed improvement we need a new algorithm.

New algorithm

One key observation is to note that once we have checked \$0145\$ we do not need to check \$1450\$ or any other combination of \$(0, 1, 4, 5)\$ since they all the the same power sum

$$ P_5[(0, 1, 4, 5)] = 0^5 + 1^5 + 4^5 + 5^5 = 4150 $$

To check whether if some permutation of \$(0, 1, 4, 5)\$ is equal to it's power sum all we need to do is to sort the digits in the power sum and compare (assuming the tuple is already sorted).

Implementation

One way to find the unique tuples is with combinations_with_replacement() from the itertools libary.

The code uses two checks to further decrease the number of combinations to check. Naively with combinations_with_replacements() there are a total of 5005 numbers to check. The check

if power_sum % 10 in perm:

is pretty intuitively. It checks whether the last number in the power_sum is in the permutations. A necessity if some permutation of perm is equal to the digit power sum. This reduces the number to check to 730. The final check is

if len_power > len_perm: return False

This cuts the number needed to check down to 591. Not this is not a huge improvement, however it is a very cheap test to perform and it is more useful for higher powers. I have attached a small table below

 P   S           Sum       # 0      # 1      # 2         time
--------------------------------------------------------------
 3   4           1301      220       60       50     0.8685 ms
 4   3          19316      715      224      165     5.0320 ms
 5   5         248860     2002      730      591    10.7200 ms 
 6   1         548834     5005     1974     1603    32.0400 ms
 7   4       25679675    11440     5040     4220    77.1300 ms
 8   3      137949578    24310    12504    10812   136.8000 ms
 9   4     2066327172    48620    24380    21224   428.4000 ms
10   1     4679307774    92378    48706    42362   858.0000 ms
11   8   418030478906   167960    92512    82456  2502.0000 ms
12   0              0   293930   175452   159471  7021.0000 ms

P is the power, S is the number of solutions. # 0 is the total number of permutations to iterate over. # 1 is the remaining numbers after the first check, and # 2 is the remaining numbers after the second check.

Code

from itertools import combinations_with_replacement as CnR

def digit_power(power):
    total = -1
    for perm in CnR(xrange(10), power):
        power_sum = sum(i**power for i in perm)
        if power_sum % 10 in perm:
            power_sum_list = map(int, str(power_sum))
            if is_fift_power(power_sum_list, perm, power):
                total += power_sum
    return total


def is_fift_power(power_lst, perm, len_perm):
    len_power = len(power_lst)
    if len_power > len_perm: return False
    power_sort = [0]*(len_perm - len_power) + sorted(power_lst)
    return power_sort == list(perm)


if __name__ == '__main__':
    import timeit
    print digit_power(5)
    times = 100
    result = timeit.timeit(
        "digit_power(5)", number=times, setup="from __main__ import digit_power")
    print 1000*result/float(times), 'ms'