Efficient search of Munchausen numbers

munchausen.py

#!/usr/bin/env python
# coding: utf8
"""
efficient search of Münchausen numbers
https://en.wikipedia.org/wiki/Munchausen_number
https://oeis.org/A046253
motivated by http://www.johndcook.com/blog/2016/09/19/munchausen-numbers/
"""

__author__ = "Philippe Guglielmetti"
__credits__ = []
__license__ = "LGPL"

from itertools import combinations_with_replacement, permutations
import datetime

numerals='0123456789abcdefghijklmnopqrstuvwxyz'

def str_base(num,b):
    return ((num == 0) and numerals[0]) or (str_base(num // b, b).lstrip(numerals[0]) + numerals[num % b])

def munchausen(base=10,zero=0):
    """generates munchhausen numbers in specified base
    :param zero: value of 0^0. if =1, numbers cannot contain zeroes
    """ 
    found=set()
    nn=[(numerals[n],n**n) for n in range(base)] # list of (nth digit,n^n)
    nn=nn[1:] #remove the 0
    for l in range(1,base):
        nnn=[(d,t) for d,t in nn if t<base**(l+1)] #consider only digits such that n^n is small enough
        for c in combinations_with_replacement(nnn,l):
            (d,t)=zip(*c) # d are the digits (without possible 0s), t are the terms
            s=sum(t) #sum of all n^n
            s0=str_base(s,base)
            s1=tuple(s0)
            s2=sorted(s1) #digits of the sum, sorted
            #removes and counts zeroes in the sum's digits
            zeroes=0
            while s2[0]=='0':
                if zero: break #numbers containing a zero are not possible if 0^0=1
                zeroes+=1
                s2.pop(0)
            #check if sum corresponds to a permutation of digits:
            if s2==sorted(d): #yes
                d=list(d)+['0']*zeroes
                for n in permutations(d):
                    if n==s1: #found !
                        if s not in found:
                            yield s0
                            found.add(s)
                            break


for base in range(2,16+1):
    t=datetime.datetime.now()
    print('base',base,':',list(munchausen(base)),datetime.datetime.now()-t)

since sum of the digits powers does not depend on their order, the algorithm consists in generating all possible combinations of the digits, (without zeroes), and check if they match the digits of the sum. If so, permutations of the digits (+ missing zeroes) are generated until we find a number corresponding to the sum.
Zeroes are handled separately in order to : 1- avoid leading zeroes, 2-handle both 0^0 conventions
Results are pretty quick until base 13, then notably slower. There's still room for improvement as shown by the "found" set which has been added because each number is found several times for an unspoted reason

results (with 0^0=0. remove numbers containing a 0 for the 0^0=1 convention):

base 2 : ['1'] 0:00:00
base 3 : ['1', '12', '22'] 0:00:00
base 4 : ['1', '130', '131', '313'] 0:00:00.001000
base 5 : ['1', '103', '2024'] 0:00:00
base 6 : ['1', '22352', '23452'] 0:00:00.001000
base 7 : ['1', '13454'] 0:00:00.004001
base 8 : ['1', '401'] 0:00:00.016001
base 9 : ['1', '30', '31', '156262', '1647063', '1656547', '34664084'] 0:00:00.0 70007
base 10 : ['1', '3435', '438579088'] 0:00:00.293030
base 11 : ['1', '66501', '517503', '18453278', '18453487'] 0:00:01.211121
base 12 : ['1', '3a67a54832'] 0:00:05.128513
base 13 : ['1', '33660', '33661', '2aa834668a', '4ca92a233518', '4ca92a233538'] 0:01:13.514350
base 14 : ['1', '23'] 0:01:29.025902
base 15 : ['1', '4b1648420dcd2', '5acbc41c19e402', '1243eed3419110e', '5acbc41c19e333', '5a99e538339a43'] 4:11:52.097591
base 16 : ['1', 'c4ef722b7820c2f', 'c76712ffc311e6e', '123b2309ff572f6b', 'c4ef722b782c26f'] 4 days, 19:39:16.014207

note the HUGE step in computation time between base 14 and base 15. I have no explanation for this ...