# Better code to find Narcissistic number

My code finding Narcissistic numbers is not that slow, but it's not in functional style and lacks flexibility: if $n \neq 7$, I have to rewrite my code. Could you give some good advice?

nar = Compile[{$}, Do[ With[{ n = 1000000 a + 100000 b + 10000 c + 1000 d + 100 e + 10 f + g, n2 = a^7 + b^7 + c^7 + d^7 + e^7 + f^7 + g^7}, If[n == n2, Sow@n]; ], {a, 9}, {b, 0, 9}, {c, 0, 9}, {d, 0, 9}, {e, 0, 9}, {f, 0, 9}, {g, 0, 9}], RuntimeOptions -> "Speed", CompilationTarget -> "C" ]; Reap[nar@0][[2, 1]] // AbsoluteTiming (*{0.398023, {1741725, 4210818, 9800817, 9926315}}*)  • 4679307774 is the sole such number of 10 digits. I'd not refer to the code implicated in this case as elegant, so much as disappointing (was hoping it could handle such a computation in less than the 15 hours it took). Oct 8, 2013 at 15:30 ## 8 Answers Here is a functional approach: Narciss[x_] := With[{num = IntegerDigits[x]}, Total[num^Length[num]] == x]  Here is a compiled version of the above function: NarcissC = Compile[{{x, _Integer}}, With[{num = IntegerDigits[x]}, Total[num^Length[num]] == x], Parallelization -> True, CompilationTarget -> "C", RuntimeAttributes -> Listable, RuntimeOptions -> "Speed"]  Now you can do something like AbsoluteTiming[Position[NarcissC[Range[10000000]], True] // Flatten]  {1.003214, {1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, 54748, 92727, 93084, 548834, 1741725, 4210818, 9800817, 9926315}} to get all the m-Narcissistic numbers from 1 to 10000000. For a further bump in speed as suggested by chyaong, here is NarcissC2 (Using Sum instead of Total) NarcissC2 = Compile[{{x, _Integer}}, With[{num = IntegerDigits[x]}, Sum[i^Length@num, {i, num}] - x], CompilationTarget -> "C", RuntimeAttributes -> Listable, RuntimeOptions-> "Speed"];  Now you can do: Pick[#, NarcissC2[#], 0] &@Range[10000000] // AbsoluteTiming  Which gives: {0.475276, {1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, 54748, 92727, 93084, 548834, 1741725, 4210818, 9800817, 9926315}} EDIT It turns out that you can get a bump using Total and Pick instead of Position (not as fast as Sum):  NarcissC1 = Compile[{{x, _Integer}}, With[{num = IntegerDigits[x]}, Total[num^Length[num]] - x], Parallelization -> True, CompilationTarget -> "C", RuntimeAttributes -> Listable, RuntimeOptions -> "Speed"]  Then Pick[#, NarcissC1[#], 0] &@Range[10000000] // AbsoluteTiming  gives: {0.626322, {1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, 54748, 92727, 93084, 548834, 1741725, 4210818, 9800817, 9926315}}  • Nice! However, your final timing uses NarcissC2 whereas you only define NarcissC Oct 7, 2013 at 8:34 • @WalkingRandomly, I was just editing that. NarcissC2 on my system is what I copied here as NarcissC Oct 7, 2013 at 8:36 • no worries. On my machine, your code takes about 2 seconds..so several times slower than the original code. How fast did the original take on yours? Ahh, I see why. You used Timing instead of Absolute Timing Oct 7, 2013 at 8:37 • You need to use AbsoluteTiming in your code for a proper comparison. Your code does more work though! Anyway, very nice. I gave it a +1 Oct 7, 2013 at 8:41 • @WalkingRandomly, you're right. Fixed it. It takes about 1 second now. so the original is about 4 to 5 times faster. But like you said, mine is doing more work. Thanks, +1 to you too. Oct 7, 2013 at 8:45 Not an answer per se, but two clarifications (which are too long for the comment box): 1) The Wiki definition you have linked to for a narcissistic number is not really apt. The Wiki page is actually describing the definition for an Armstrong Number, also known as pluperfect digital invariants, or m-narcissistic numbers, such as: $$407 = 4^3 + 0^3 + 7^3$$ These require the use of the power term$m$(the 3 in this example) over and above the digits of integer$n= 407$. By contrast, the correct and proper reference to the term 'narcissistic number' comes from the article by Madachy, J. S. (1966), Mathematics on Vacation, Thomas Nelson & Sons — p.163 to 175, who defines them as numbers: "that are representable, in some way, by mathematically manipulating the digits of the numbers themselves". What the Wiki page describes ... the Armstrong numbers ... is quite different, ... not the 'narcissistic numbers', as the page claims, but the m-narcissistic numbers. But that's wiki for you. 2) Finite to infinite: The set of narcissistic numbers involve a finite search (such as the solutions above). The problem becomes rather more tricky if you allow for the use of radicals or factorials ... because the search problem is no longer finite ... rather you can have infinite nesting of square root symbols or factorial symbols. One can get some pretty results when you allow radicals, such as say: For more detail, please see: http://www.tri.org.au/numQ/pwn/ or a fun little piece I did entitled: Radical Narcissistic Numbers, Journal of Recreational Mathematics, 33(4), 2004-2005, 250-254. I've been meaning to put up the mma code for this too ... this was done long before the age of multi-processors, so I think I'll have to update the code for parallel cores, which would make an enormous difference here. nar[m_] := ToExpression[ "Compile[{$},Do[With[{n=0" <>
StringJoin[
Table["+1" <> Array["0" &, m - 1 - i, 1, StringJoin] <> "a" <>
ToString[m - 1 - i], {i, 0, m - 1}]] <> ",n2=0" <>
Table["+a" <> ToString[m - 1 - i] <> "^" <> ToString[m], {i, 0,
m - 1}] <> "},If[n\[Equal]n2,Sow@n];];,{a0" <>
StringJoin[Table[",0,9},{a" <> ToString[i], {i, 1, m - 1}]] <>
",9}],RuntimeOptions\[Rule]\"Speed\",CompilationTarget\[Rule]\"C\"\
]"];

Reap[nar[7][0]][[2, 1]] // AbsoluteTiming

(*{1.184733, {9926315, 1741725, 9800817, 4210818}}*)


My computer is rather slow. @RunnyKine's code takes 0.901549 seconds on my computer.

• Interesting way of code generation :) Oct 7, 2013 at 15:38
• @IstvánZachar I'm thinking the same!
– CHM
Oct 15, 2013 at 15:42

From a cold start, I would have written it like this:

 findNarc = Compile[{{stop, _Integer}, {pow, _Integer}},
Do[
If[Total[IntegerDigits[n]^pow] == n, Sow[n]]
, {n, 1, stop}
]
, RuntimeOptions -> "Speed", CompilationTarget -> "C"];


However, it is slower than your function (which takes 0.326 seconds on my machine)

Reap[findNarc[10000000, 7]] // AbsoluteTiming

(*{2.900166, {Null, {{1, 1741725, 4210818, 9800817, 9926315}}}}*)


I'd usually use InternalBag instead of Sow since it can be compiled whereas Sow cannot (see InternalBag inside Compile ) but Sow is called so infrequently here that I don't think that's the problem. Besides, you used it as well so both my code and your code would have been hit by the same penalty.

So, for the sake of speed, I'd be tempted to just do a very minor modification of your code to give the flexibility for choosing the Power:

nar =
Compile[{{pow, _Integer}},
Do[With[{n =
1000000 a + 100000 b + 10000 c + 1000 d + 100 e + 10 f + g,
n2 = a^pow + b^pow + c^pow + d^pow + e^pow + f^pow + g^pow},
If[n == n2, Sow@n];];, {a, 9}, {b, 0, 9}, {c, 0, 9}, {d, 0,
9}, {e, 0, 9}, {f, 0, 9}, {g, 0, 9}], RuntimeOptions -> "Speed",
CompilationTarget -> "C"];

Reap[nar[7]] // AbsoluteTiming

(*{0.329019, {Null, {{1741725, 4210818, 9800817, 9926315}}}}*)


Dynamically generated Do loops:)

cnar =
With[{n = 7},
With[{var = Array[Unique["x"] &, n]},
With[{n1 = FromDigits@var, n2 = Total[var^n]},
Compile[{Null},
Do[If[n1 == n2, Sow@n1], ##],
RuntimeOptions -> "Speed", CompilationTarget -> "C"
] & @@ MapAt[1 &, Thread[{var, 0, 9}], {1, 2}]
]
]
];

Reap[cnar@0][[2, 1]] // AbsoluteTiming

(*
CompiledFunction[{},
Do[If[10 (10 (10 (10 (10 (10 x3 + x4) + x5) + x6) + x7) + x8) + x9 ==
x3^7 + x4^7 + x5^7 + x6^7 + x7^7 + x8^7 + x9^7,
Sow[10 (10 (10 (10 (10 (10 x3 + x4) + x5) + x6) + x7) + x8) +
x9]], {x3, 1, 9}, {x4, 0, 9}, {x5, 0, 9}, {x6, 0, 9}, {x7, 0,
9}, {x8, 0, 9}, {x9, 0, 9}], "-CompiledCode-"]

{0.358024, {1741725, 4210818, 9800817, 9926315}}
*)


It should be much faster to generate all possible integer digit sets, and then select those integer digit sets that have the require property. For instance, $135, 153, 315, 351, 513, 531$ all have the integer digits $1, 3, 5$, but the sum of the cubes of all 6 digit sets is the same, namely, $153$. The set of all possible integer digit sets for an integer length $k$ number (assuming leading 0s are acceptable) is $\binom{k+9}{k}$, which is much smaller than the set of all integers with integer length less than or equal to $k$. For instance, for $k=7$, we have:

10^7-1
Binomial[9+7, 7]


9999999

11440

Checking each of the integer digit sets will be much faster than checking all integers. Here is a function that carries out this idea:

narcisist[p_, k_] := With[
{subs=Transpose[Transpose @ Subsets[Range[0, 8+k], {k}] - Range[0, k-1]]},

With[
{tdigits=Sort/@IntegerDigits[Total[subs^p,{2}],10,k]},

Total[Pick[subs,Total[Abs[subs-tdigits],{2}],0]^p,{2}]
]
]


I extended the concept to having independent powers $p$ and integer lengths $k$. Here is the answer for the OP parameters:

narcisist[7, 7] //AbsoluteTiming


{0.020188, {0, 1, 9800817, 4210818, 1741725, 9926315}}

This approach can succeed for much larger integer lengths:

narcisist[10, 10] //AbsoluteTiming


{0.215505, {0, 1, 4679307774}}

narcisist[20, 20] // AbsoluteTiming


{112.554, {0, 1, 63105425988599693916}}

A fast method：

nar[n_] := Pick[#, #~BitXor~Range[10^(n - 1), 10^n - 1], 0] &@
Flatten[Outer[Plus, ##] & @@ Array[Range[Boole[# == 1], 9]^n &, n]]

nar /@ Range[7] // AbsoluteTiming


{0.314020, {{1, 2, 3, 4, 5, 6, 7, 8, 9}, {}, {153, 370, 371, 407}, {1634, 8208, 9474}, {54748, 92727, 93084}, {548834}, {1741725, 4210818, 9800817, 9926315}}}

This may not be efficient but it is terse:

narc[n_] := Module[{r, l, t},
r = Range[n];
t = Total@(#^Length[#]) & /@ (IntegerDigits /@ r);
Pick[r, r - t, 0]
]


narc[10000000] yields:

{1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, \
54748, 92727, 93084, 548834, 1741725, 4210818, 9800817, 9926315}

• A bit faster：Module[{r, l, t}, r = Range[10^6]; t = Total[With[{id = IntegerDigits@#}, id^Length@id]] & /@ r; Pick[r, r~BitXor~t, 0]] // AbsoluteTiming` Oct 8, 2013 at 12:14
• @chyaong thank you...nice Oct 8, 2013 at 12:17