How to Generate Münchhausen Numbers in High Radices? A Münchhausen number is a whole number equal to the sum of its digits raised to powers of themselves. For the purpose of such calculations, the convention is that $0^0 = 1$. For example, in radix 10, the number 3435 is a Münchhausen number, since $3435 = 3^3 + 4^4 + 3^3 + 5^5$. In radix 4, the number 313 is a Münchhausen number, since $313 = 3^3 + 1^1 + 3^3$. On a curious property of 3435 covers more details, and shows that Münchhausen numbers in radix $b$ must be less than or equal to $2b^b$.

Münchhausen numbers for base $b$ can be generated with a brute-force search, by checking whether each number from 1 to $2b^b$ is equal to the sum of its digits raised to powers of themselves. One way to improve upon this is to generate multisets of digits, and for each combination check whether the sum of digits raised to powers of themselves produces a number with the same multiset of digits. Even with parallelization, this seems infeasible for large bases. For example, for base 35, a Münchhausen number would have at most 36 digits. The number of combinations (with repetitions) of 36 elements from 35 digits is $\binom{35 + 36 - 1}{36}=\binom{70}{36}=109,069,992,321,755,544,170.$ This is in addition to those multisets with fewer than 36 elements that would also be checked.

If any, what algorithm(s) can be used to generate Münchhausen numbers feasibly for high bases (e.g., base 35)? OEIS sequence A166623 includes Münchhausen numbers for bases between 2 and 35.
Meet-in-the-middle

One approach is to use a meet-in-the-middle algorithm.  This will reduce the running time from  $2b^b$ to something in the ballpark of $\sqrt{2b^b}$.

In particular, let your candidate number be
$$A = a_b b^b + \dots + a_1 b + a_0.$$
We'll write this in the form $A=A^* + A^\dagger$ where
$$\begin{align*}
A^* &amp;= a_b b^b + \dots + a_{\lfloor b/2 \rfloor+1} b^{\lfloor b/2 \rfloor+1}\\
A^\dagger &amp;=a_{\lfloor b/2 \rfloor} b^{\lfloor b/2 \rfloor} + \dots + a_0.
\end{align*}$$
Finally, let $s(a_i b^i + \dots + a_0) = a_i^{a_i} + \dots + a_0^{a_0}$.  We'll look for all solutions to the equation
$$A^* + A^\dagger = s(A^*) + s(A^\dagger).$$
This is equivalent to
$$A^* - s(A^*) = s(A^\dagger) - A^\dagger.$$
We'll compute a table of all possible values for $s(A^\dagger)-A^\dagger$, where $A^\dagger$ varies over all possibilities of the form above.  Then, we'll compute $A^*-s(A^*)$, for each $A^*$ of the form above, and look it up in the table to check for a match.  Each match yields a solution to the equation above, and thus a valid Münchhausen number.  This can be used to enumerate all Münchhausens.

For $b=10$, the running time becomes something like $10^5$ steps of computation, which should be no problem.  However, this is entirely infeasible for $b=35$.  The algorithm is still exponential-time.

Linear programming

Another possible idea is to use integer linear programming.

Let $x_{i,d}$ be a 0-or-1 variable, with the intended meaning that $x_{i,d}=1$ means that $a_i = d$.  Then we can obtain a bunch of linear inequalities that characterize a correct solution.  In particular, we have
$$\sum_{i,d} x_{i,d} d b^i = \sum_{i,d} x_{i,d} d^d.$$
We also have $\sum_d x_{i,d} = 1$ for each $i$.

Next, given a lower and upper bound $L,U$, we could use integer linear programming to search for a valid solution in the range $[L,U]$, by adding the inequality
$L \le \sum_{i,d} x_{i,d} d b^i \le U$ and looking for a valid solution.  To enumerate all solutions, start with a broad range $[0,2b^b]$, then when you find a valid solution, split the current range in two and recurse.  Then, use an off-the-shelf ILP solver to find solutions.

I don't expect this to be very effective, but you could augment it with additional inequalities that might help the ILP solver run faster.  In particular, for each prime number $p$ less than $2b^b$, we have
$$\sum_{i,d} x_{i,d} d b^i \equiv \sum_{i,d} x_{i,d} d^d \pmod p$$
This can be encoded via the linear equality
$$\sum_{i,d} x_{i,d} d b^i = \sum_{i,d} x_{i,d} x_{i,d} d^d + y_p p$$
where $y_p$ is an integer variable.  Note that $d b^i$ and $d^d$ are integer constants, so you can reduce them modulo $p$. Thus we get the linear equality
$$\sum_{i,d} \alpha_{i,d,p} x_{i,d} = \sum_{i,d} \beta_{d,p} x_{i,d} + y_p p$$
where $\alpha_{i,d,p} = d b^i \bmod p$ and $\beta_{d,p} = d^d \bmod p$ are constants you can precompute.  You can do the same if $p$ is a power of a prime.

So, you might add linear inequalities of this form for several well-chosen prime powers (say, all of the prime powers below some bound; or all small prime numbers, plus the prime power factors of $b$, $b-1$, $b+1$, $b^2$, $b^2-1$, $b^2+1$).  One could hope that these help the linear solver find solutions faster.

I don't know whether this will be effective, but it is something you could try.  It would not be too hard to implement and see how well it works.