In today’s Programming Praxis exercise, our task is to calculate all the narcissistic numbers, also known als the Armstrong numbers or the pluperfect digital invariants, i.e. the sequence of numbers for which the sum of the cubes of the digits is equal to the number itself.

Supposedly, a mathematician by the name of Dik Winters developed an algorithm in 1985 that could generate all 88 numbers in about half an hour, which should theoretically run in seconds on modern day hardware. Unfortunately, neither Phil (the author of the Programming Praxis blog) nor I were able to find the original algorithm. The exercise provides a brute-force solution, which of course will not terminate in anything close to an acceptable time since the highest number in the sequence has 39 digits.

The solution I came up with in the end (after lots of different approaches to speed things up) is significantly faster than the naive brute-force solution, yet still nowhere close to the theoretical solution.

Some imports:

import Data.List
import qualified Data.Vector as V

Since calculating the full sequence takes too long, I added an argument to specify the maximum amount of desired digits for timing purposes.

narcissistic :: Integer -> [Integer]
narcissistic upto = narcs =<< [1..min 39 upto] \\ [2,12,13,15,18,22,26,28,30,36]

When generating the narcissistic numbers, we make a number of improvements to the brute-force algorithm:

- Since the order of the digits doesn’t matter for the sum of cubes, we only look increasing series of digits, ruling out all permutations
- Since the power function is relatively expensive, we precalculate all 10 possibilities into a lookup table
- All digits sequences with a sum that is too low or high for an n-digit number are ignored

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)

With this approach, calculating the numbers of 1 through 16 digits takes about 3.4 seconds. 1 through 25 digits takes just over 4 minutes. I have no idea how long the full sequence would take. I’m sure there’s some way to eliminate more options using some mathematical proof, but I haven’t been able to find or come up with one.

main :: IO ()
main = mapM_ print $ narcissistic 16