Here is a nice mathematical exercise from Programming Praxis. Is about finding the “narcissistic numbers”, n digit numbers numbers where the sum of all the nth power of their digits is equal to the number.

To reduce the problem a little, I decided to start by limiting the number of digits. So, the first approach will be just calculate if a number is narcissistic of not. So, after checking it and making a couple of performance adjustments, the code is as follows…

import math
NUM_DIGITS = 8

powers = []
for i in xrange(NUM_DIGITS + 1):
    powers.append([j ** i for j in xrange(10)])

def number_is_narcissistic(number):
    str_number = str(number)
    num_digits = len(str_number)
    power = powers[num_digits]
    acc = 0
    for i in (int(d) for d in str_number):
        acc += power[i]
        if acc > number:
            return False

    return acc == number

def number_is_narcissistic_no_str(number):
    acc = 0
    original_number = number
    num_digits = int(math.log10(number)) + 1
    power = powers[num_digits]
    while number:
        acc += power[number % 10]
        number = number / 10
        if acc > original_number:
            return False

    return acc == original_number

def main():
    ''' Find all the narcissistic numbers'''
    for i in xrange(1, 10 ** NUM_DIGITS):
        if number_is_narcissistic_no_str(i):
            print i

if __name__ == '__main__':
    main()

Some comments:

• To get the digits of a number, the more Pythonic way is to convert the number to string and then get each of the chars (converting them to integers). After testing, that’s not as fast as using a loop and divide a number by 10 several times. I left the original function, though.
• The same thing happens to calculate the number of digits. Again, I prefer to use len(str(number)) as is easier to read and understand, but calculating the log is faster.
• Also, precalculating the powers of each digit an storing them (in the ‘powers’ array) is faster than calculating them all the time.
• Adding a ‘return’ statement makes us exit the loop as soon as there is no possibility of making it.

With all those optimizations, the time to calculate all the narcissistic numbers up to 8 digits is around 4m 40s in my laptop.

That seems a little long, even if we are in the Python performance realm. So, I scratched my head a little more and get with this idea. Any possible combination of digits, no matter how is ordered, can only produce one result. For example, the combination 153 (which is a narcissistic number) has the same digits than 531 or 315. But, if we calculate the result, we can check only the result as a possible narcissistic number. So, with that idea in mind, I made the program again, this time getting all the possible combinations for  number of digits and obtaining a list of candidates that is way lower than ever possible number.

from itertools import combinations_with_replacement
import math
NUM_DIGITS = 8

powers = []
for i in xrange(NUM_DIGITS + 1):
    powers.append([j ** i for j in xrange(10)])

def number_is_narcissistic_no_str(number):
    acc = 0
    original_number = number
    num_digits = int(math.log10(number)) + 1
    power = powers[num_digits]
    while number:
        acc += power[number % 10]
        number = number / 10
        if acc > original_number:
            return False

    return acc == original_number

def comb_is_candidate(comb, num_digits):
    ''' Combination of number is candidate, as his power has the same number of digits '''
    acc = 0
    limit = 10 ** num_digits
    min_limit = 10 ** (num_digits - 1)
    power = powers[num_digits]
    for n in comb:
        acc += power[n]
        if acc > limit:
            return False

    return acc > min_limit

def number_power(comb, power):
    return sum(c ** power for c in comb)

def all_candidate_numbers(num_digits):
    comb = combinations_with_replacement(xrange(10), num_digits)
    # Check the possible numbers
    candidate_groups = [number for number in comb if comb_is_candidate(number, num_digits)]
    candidates = set()
    for candidate in candidate_groups:
        candidates.add(number_power(candidate, num_digits))

    return candidates

def main():
    for digits in xrange(1, NUM_DIGITS + 1):
        candidates = all_candidate_numbers(digits)
        for candidate in candidates:
            if number_is_narcissistic_no_str(candidate):
                print candidate

if __name__ == '__main__':
    main()

Comments to this version:

• Using itertools module, generates all the possible digits combinations. Then, get if that combination can produce a narcissistic number. For example, for 8 digits, the number of candidates is 16,131. Quite lower than 90,000,000 numbers.
• Please note that the program is not checking if the result of a combination is itself a narcissistic number. It could have other digits. But it uses the result of powering all the digits as a candidate, and checks the result again. Some more optimisations could be investigated on that area.
• One interesting difference is that, while it will print all the narcissistic numbers, as the first version does, it won’t print them in the same order, as the generation of the candidates is not generating them from lower to higher. If that’s an issue, the results or the candidates can be sorted.
• The rest of the comments of the first version applies (not converting to string, storing the powers, etc)

So, that gives a time of … 450 ms to get up to 8 digits combinations! That’s quite an improvement. I pushed it a little and get  about 3m 40s for getting all the numbers up to 18 digits… The highest possible narcissistic number has 39 digits, but I’m not going to get that far… 😉 I think that’s quite good result for Python, and I am quite happy with the result.