115132219018763992565095597973971522401
December 14, 2012
The digits and sum functions of the Standard Prelude make this easy:
(define (pow n) (lambda (x) (expt x n)))
> (do ((n 1 (+ n 1))) (#f)
(let ((ds (digits n)))
(when (= (sum (map (pow (length ds)) ds)) n)
(display n) (newline))))
1
2
3
4
5
6
7
8
9
153
370
371
407
1634
8208
9474
54748
92727
93084
548834
1741725
4210818
9800817
9926315
24678050
24678051
88593477
146511208
472335975
...
You can run the program at http://programmingpraxis.codepad.org/A1240SS8.
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?
Programming Praxis 
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.
import Data.List import qualified Data.Vector as V narcissistic :: Integer -> [Integer] narcissistic upto = narcs =<< [1..min 39 upto] \\ [2,12,13,15,18,22,26,28,30,36] narcs :: Integer -> [Integer] narcs n = sort $ f [] 0 9 n where powers = V.fromList $ map (^n) [0..9] pow i = powers V.! fromIntegral i (lo, hi) = (10^(n-1), 10^n) f ds s x 1 = [ s' | i <- [0..x], let s' = s + pow i, s' >= lo , s' < hi, sort (show s') == (show =<< (i:ds))] f ds s x d = [0..x] >>= \i -> f (i:ds) (s + pow i) i (d-1)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:
def find_narc(limit): n=1 i=[] while n<=limit: k=0 number=str(n) l=len(number) for x in range(0,l): k=k+int(number[x])**l if n==k: i.append(n) n=n+1 return i [/sourcecode ] Which gives the following results for a limit 10,000,000: [sourcecode lang="python"] print find_narc(10000000) # [1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, 54748, 92727, 93084, 548834, 1741725, 4210818, 9800817, 9926315] # total time 67.6sI 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:
def find_narc2(limit): ret=[] i=[0,1,2,3,4,5,6,7,8,9] for z in range(1,limit+1): for x in range(0,10): i[x]=x**z for n in range(10**(z-1),(10**z-1)): number=str(n) k=0 for y in range (0,z): k=k+i[int(number[y])] if n==k: ret.append(n) return retWhen 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):
from itertools import permutations,combinations_with_replacement # Returns a list of all unique N-digit combinations (single permutation) def uniques(N): return [''.join(l) for l in list(combinations_with_replacement('0123456789',N))] # Returns a list of all number permutations of the digit-string n def permut(n): pm = permutations(str(n)) return [int(''.join(tpl)) for tpl in set(pm)] def fastnarcs(N): narcs = [] for k in range(1,N+1): for n in uniques(k): S = sum([int(i)**k for i in list(n)]) for pmt in permut(n): if int(pmt) == S: #print '%s gave narc %d'%(n, pmt) narcs.append(pmt) narcs = list(set(narcs)) narcs.sort() return narcs def simplenarcs(N): narcs = [] for n in range(0,10**N): k = len(str(n)) S = sum([int(i)**k for i in list(str(n))]) if S == n: narcs.append(n) return narcsfastnarc 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 […]
If you are interested in number theory and recreational mathematics, and wish to practice your programming skills further, indulge yourself with these number curiosities: http://www.glennwestmore.com.au/category/advanced-number-curiosities/ or http://www.glennwestmore.com.au/category/number-curiosities/
I revisited this problem recently and have relatively simple Java code to solve it in ~1.2 seconds on my laptop (searching all lengths up to 60), plus a base-16 version that takes around 60 hours (searching all lengths up to 116): https://github.com/peterdettman/armstrong-numbers .