9
$\begingroup$

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}}*)
Improve this question
$\endgroup$
1
  • $\begingroup$ 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). $\endgroup$
    – Daniel Lichtblau
    Oct 8, 2013 at 15:30

8 Answers 8

Reset to default
10
$\begingroup$

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}}
Improve this answer
$\endgroup$
9
  • $\begingroup$ Nice! However, your final timing uses NarcissC2 whereas you only define NarcissC $\endgroup$
    – WalkingRandomly
    Oct 7, 2013 at 8:34
  • $\begingroup$ @WalkingRandomly, I was just editing that. NarcissC2 on my system is what I copied here as NarcissC $\endgroup$
    – RunnyKine
    Oct 7, 2013 at 8:36
  • $\begingroup$ 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 $\endgroup$
    – WalkingRandomly
    Oct 7, 2013 at 8:37
  • $\begingroup$ 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 $\endgroup$
    – WalkingRandomly
    Oct 7, 2013 at 8:41
  • $\begingroup$ @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. $\endgroup$
    – RunnyKine
    Oct 7, 2013 at 8:45
10
$\begingroup$

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.

enter image description here

One can get some pretty results when you allow radicals, such as say:

enter image description here

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.

Improve this answer
$\endgroup$
6
$\begingroup$ 
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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Interesting way of code generation :) $\endgroup$
    – István Zachar
    Oct 7, 2013 at 15:38
  • $\begingroup$ @IstvánZachar I'm thinking the same! $\endgroup$
    – CHM
    Oct 15, 2013 at 15:42
5
$\begingroup$

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 Internal`Bag 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}}}}*)
Improve this answer
$\endgroup$
5
$\begingroup$

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}}
*)
Improve this answer
$\endgroup$
4
$\begingroup$

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

Improve this answer
$\endgroup$
3
$\begingroup$

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

Improve this answer
$\endgroup$
2
$\begingroup$

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}
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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 $\endgroup$
    – chyanog
    Oct 8, 2013 at 12:14
  • $\begingroup$ @chyaong thank you...nice $\endgroup$
    – ubpdqn
    Oct 8, 2013 at 12:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.