# 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). – Daniel Lichtblau Oct 8 '13 at 15:30 ## 7 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 – WalkingRandomly Oct 7 '13 at 8:34 @WalkingRandomly, I was just editing that. NarcissC2 on my system is what I copied here as NarcissC – RunnyKine Oct 7 '13 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 – WalkingRandomly Oct 7 '13 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 – WalkingRandomly Oct 7 '13 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. – RunnyKine Oct 7 '13 at 8:45 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}}}}*)  - 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 :) –  István Zachar Oct 7 '13 at 15:38
@IstvánZachar I'm thinking the same! –  CHM Oct 15 '13 at 15:42

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}} *)  - 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:

http://www.numq.com/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.

-

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 –  chyaong Oct 8 '13 at 12:14
@chyaong thank you...nice –  ubpdqn Oct 8 '13 at 12:17

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}}}

-