The question referred to at the beginning of the exercise was how to write a program that calculates narcissistic numbers. Clearly there can be no narcissistic numbers greater than sixty digits, since n · 9n < 10n−1 for n > 60. Wolfram MathWorld mentions that Dik Winter calculated the complete list of narcissistic numbers in 1985; I looked up the reference, and found that it took 2000 seconds (about half an hour) on a Vax, which probably equates to a second or two on a recent-vintage PC, but I don’t know how to make that calculation; certainly indexing through the integers as done above will take far too long. Can anybody help?

Well, it certainly doesn’t run in a few seconds but it’s at least significantly faster than the direct solution. It’s calculated all of the Narcissistic numbers up to the five 25 digit solutions in about 45 minutes and it’s still running. It’s took about 30 seconds to do what Mithrandir’s similar algorithm did in 1 though, so there’s room for improvement.

I’ll see if I can think of any more optimizations ( or implement them if anyone else has some :) ), but that’s what I have for now.

Just a thought: Let us say you have tested all permutations of N digits. Now you add a digit 0..9 number N+1 to the previous digits and calculate a new number. Does this calculated number have a permutation of digits equal to the digits including the new one? If it does, the calculated number tells the correct permutation of digits. Try next digit. If no digits 0..9 gives a new narcissistic number then ?

I feel like there is a better way to do this geometrically, since you’re computing sums of digits to the nth power. Maybe a predictive solution where the computer computes partial sums and then finds integer n-roots for remaining difference.

It generates the possible candidates checking all the digits combinations using itertools, then checking again those candidates to get the narcissistic numbers. It seems quite fast to me (for being in Python), around 3m 40s for all numbers up to 18 digits.

I have written a solution to this problem (in C#, available as a VS 2012 Express solution here). The code is somewhat involved, but it searches the entire problem space and finds the 88 solutions other than “0”. Running time is 1.6 seconds as a single-threaded program on an i7-2640M @2.80GHz. So I believe I can now answer the question posed by the post.

The code is in 3 parts: a base-10 multi-precision implementation (Base10e8Nat), the main Solver class that tackles each length independently, and a top-level Program that collects counts and timings for each length. A few tables are pre-calculated, among them the powers themselves and all the possible multiples of them. There is no multiplication during the actual search. Counts of powers are packed into a 64-bit integer, 6 bits for each digit count.

The search is done on the multiples of powers rather than the integers, as others have suggested. However, it begins with the 9s and recurses to 8s and so on down. At each level, we begin with the largest multiple that will not overflow the result, and then work our way down to 0, recursing for each. At each recursive search, we have available the number of each (higher) power already chosen (and the number remaining available to us), and a running total of the sum of those powers so far.

The virtue of searching this way (besides reducing it to a combinatoric search) is that as we reach the lower levels (about 6 and below for the larger lengths) of the recursion, some (or many) of the most-significant digits of the running total will already be “fixed” (modulo a single overflow), because we have a hard bound on the number of least-significant digits that could potentially change during this recursion (it is the number of digits in the highest multiple of this level’s power that are still available). We can now look at those fixed digits (hereafter the “prefix”) and apply heuristics to see if they are consistent with the already-chosen powers, and the remaining choices available.

There are two main heuristics, applied with a slight variation to the higher powers and the lower powers. Essentially there’s either too many of a digit in the prefix, or not enough.These heuristics must also take note of the possibility of overflow, but it’s not as bad as one might think at first. It has potentially unbounded effects for 9s and 0s, but the counts of other digits can only be affected by at most 1 by any overflow. We look at what the overflow digit (the LSD of the fixed part) is and in a pinch we make a conservative assumption to avoid analyzing further.

For the higher powers: the number of each power chosen by higher levels (the multiple) is fixed as far as this sub-search is concerned. So if there are already more than that number of the digit in the prefix, it can’t be fixed and this sub-search is a bust. Likewise if the digit count is lower than the multiple for some N, it will have to be made up in the remaining digits (the suffix). Observe though that if there are undercounts for several Ns, there have to be enough remaining digits to cover _all_ of them. So we can sum up the “excess powers” and then compare it to the length of the suffix. If the suffix is too short to supply enough of the missing digits, it’s a bust.

For the lower powers: similar to “excess higher powers”, here we haven’t yet apportioned the number of powers, but we do have a limit on how many are available for this sub-search. So we look at the prefix digits and figure out how many of the lower powers are already spoken for. If it’s more than we’ve got, then bust. The converse, I haven’t implemented yet: if there are only a few digits of lower powers, we shouldn’t bother recursing with too many “available” powers.

0, 9, and the current level require somewhat extra attention, which I won’t go into: performance becomes much more tractable with just a few of these heuristics applied anyway, though it really pays off to wind them nice and tight. To give you an idea of the effect, from profiling results I can see that this search only makes a full comparison of digit counts with power counts for about 1.14 million candidates. Everything else is pruned.

That’s probably enough to get the general approach across. I’ve tried to add explanatory comments through the interesting parts of the code if anyone’s curious, or ask questions here and I will try to answer.

I forgot to add: if you want to watch the recursive search in action, you can remove the Conditional attribute from Solver.DumpPowerBits(). Each of the 1.14 million “final candidates” will then be printed out.