## 115132219018763992565095597973971522401

### December 14, 2012

Today we have an exercise, a solution, and a question I don’t know the answer to. The exercise is to write a program that lists the *narcissistic* numbers in which the sum of the *n*th powers of the digits of an *n*-digit number is equal to the number. For instance, 153 is a narcissistic number of length 3 because 1^3 + 5^3 + 3^3 = 1 + 125 + 27 = 153. The largest decimal narcissistic number is the 39-digit number that is the title of today’s exercise.

Your task is to write a program to find narcissistic numbers. When you are finished, you are welcome to read or run a suggested solution, or to post your solution or discuss the exercise in the comments below.

Pages: 1 2

Partially better solution is to generate non-decreasing sequences of digits and compute their sum. If this sum has the same digits as the digits used to generate it, you have a number.

A quickly implemented Haskell solution follows. It could be optimized a lot. On my machine, it prints up to 4338281769391371 in 1 second. I don’t have time to profile it though.

import Data.List

digits :: [Int]

digits = [0..9]

expand :: Int -> [[Int]] -> [[Int]]

expand 0 l = l

expand n [] = expand (n – 1) $ map (:[]) digits

expand n l = expand (n – 1) $ concatMap exp l

where

exp ll@(x:xs) = map (:ll) $ filter (>= x) digits

test :: Int -> [Int]

test n = map (toSum 0 n) . filter (narcis n) . expand n $ []

narcis :: Int -> [Int] -> Bool

narcis n l = let x = toSum 0 n l in sort l == sort (toList x)

toSum :: Int -> Int -> [Int] -> Int

toSum s _ [] = s

toSum s e (x:xs) = toSum (s + x ^ e) e xs

toList :: Int -> [Int]

toList 0 = []

toList n = (n `mod` 10) : toList (n `div` 10)

printOne n = mapM print $ test n

main = do

mapM printOne [1..]

(hopefully WordPress will keep the code as it is)

(let’s try again, please delete my previous comment and strip out this line:P)

Partially better solution is to generate non-decreasing sequences of digits and compute their sum. If this sum has the same digits as the digits used to generate it, you have a number.

A quickly implemented Haskell solution follows. It could be optimized a lot. On my machine, it prints up to 4338281769391371 in 1 second. I don’t have time to profile it though.

import Data.List

digits :: [Int]

digits = [0..9]

expand :: Int -> [[Int]] -> [[Int]]

expand 0 l = l

expand n [] = expand (n - 1) $ map (:[]) digits

expand n l = expand (n - 1) $ concatMap exp l

where

exp ll@(x:xs) = map (:ll) $ filter (>= x) digits

test :: Int -> [Int]

test n = map (toSum 0 n) . filter (narcis n) . expand n $ []

narcis :: Int -> [Int] -> Bool

narcis n l = let x = toSum 0 n l in sort l == sort (toList x)

toSum :: Int -> Int -> [Int] -> Int

toSum s _ [] = s

toSum s e (x:xs) = toSum (s + x ^ e) e xs

toList :: Int -> [Int]

toList 0 = []

toList n = (n `mod` 10) : toList (n `div` 10)

printOne n = mapM print $ test n

`main = do`

mapM printOne [1..]

(hopefully WordPress will keep the code as it is)

I give up. Don’t know how to format it properly :|

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.

Blog post:

Narcissistic Numbers

Source code:

narcissistic.rkt

[...] today’s Programming Praxis exercise, our task is to calculate all the narcissistic numbers, also known als [...]

Here’s my Haskell solution (see http://bonsaicode.wordpress.com/2012/12/14/programming-praxis-115132219018763992565095597973971522401/ for a version with comments). Sadly it’s not the mythical 2-second solution, but it’s reasonably quick. Calculating all numbers of 1 through 25 digits takes just over 4 minutes.

My first version in Python. Still a little slow. 17 minutes for n <= 25.

My java solution here: hamidj.wordpress.com

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 ?

[...] Pages: 1 2 [...]

Here’s a simple python script to inspect all numbers to a defined limit:

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.

Well I screwed that sourcecode up.

Rather than take each digit to the nth degree for every single number, I tried to break it out by power series:

When compared to the function above, it turns out it takes about a minute (8.6%) off of 100,000,000 numbers:

http://www.solveet.com/exercises/Emulando-a-Dik-T-Winter/155/solution-1073

N=20 ___ 1,91 secs

N=25 ___ 13,12 secs

N=39 ___ 13,06 minutes

Nice solution! Loop unrolling works great here.

Simple and slightly faster variant in python (cheating with itertools):

fastnarc is about 4 times as fast as the trivial solution on N<=7

[...] is a nice mathematical exercise from Programming Praxis. Is about finding the “narcissistic numbers”, n digit numbers numbers where the sum of [...]

I’ve written a solution on my blog (in Python) http://wrongsideofmemphis.com/2012/12/18/narcissistic-numbers/

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.

[...] – As mentioned in the original post Dik Winter calculated all narcissistic numbers in 1985, which took 30 minutes then (equivalent to [...]