Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

In recreational number theory, a narcissistic number (also known as a pluperfect digital invariant (PPDI), an Armstrong number(after Michael F. Armstrong) or a plus perfect number) is a number that is the sum of its own digits each raised to the power of the number of digits.

For example,the three digits armstrong number:

enter image description here

Now I want to find the 21 digits armstrong number.I know the answer is {449177399146038697307, 128468643043731391252}.

Better code to find Narcissistic number.In this post,they use brute force method to test all number.It is impossible to use this method to find 21digits armstrong number.

My thought is just consider the count of numbers. For example, there are two 2, two 7,one 9.Since 2*2^5+2*7^5+1*9^5=92727,so {{2,2,1},{2,7,9}} is correct combination and 92727 is 5 digits armstrong number.

My code is:

Clear["Global`*"]
n = 21;
counts = IntegerPartitions[n, 10];
test1 = Compile[{{counts, _Integer, 1}, {num, _Integer, 1}}, 
    10.^(n - 1) < counts.num^(1. n) < 10.^n];
test2 = Compile[{{list, _Integer, 1}, {counts, _Integer, 1}, {num, _Integer, 1}}, 
   Catch[
     Do[If[Count[list, num[[i]]] != counts[[i]], Throw[False]], {i, Length[counts]}];     
     Throw[True]]];
f[counts_] := Module[{allpossible, allnum},
  allpossible = Permutations@counts;
  allnum = Subsets[Range[0, 9], {Length[counts]}];
  Select[Join @@ Outer[List, allpossible, allnum, 1], 
     test1 @@ # && (test2[IntegerDigits[#1.#2^n], #1, #2] & @@ #) &]
  ]
Length@counts
Monitor[#1.#2^n & @@@ Join @@ Table[f@counts[[i]], {i, %}], 1.i/%] // AbsoluteTiming

It will take about 6~10 minutes to find 21digits armstrong number.But the maximal armstrong number is 39 digits, so my code is not enough high efficiency.How can I increase of efficiency?

share|improve this question
    
Way more work than needed - cut your search space to all the non-decreasing sequences of N digits, check these for matching sums after raising to power (pre-compute the latter), if sum sorted non-decreasing matches sequence, you have a hit... –  rasher Jul 9 at 7:05
    
@rasher I'm afraid that's not sufficiently efficient. For N digits you'll need to sift through Binomial[N+9,9] non-decreasing sequences, which gives you 1.677.106.640 combinations for N = 39. –  Teake Nutma Jul 11 at 8:55
    
@TeakeNutma: Add parity check, some other rudimentary sieving, you're at a few hundred million checks for all the base 10 cases. Should run in reasonable time. In any case, +1 to the OP offering a bounty on a "look it up on OEIS" kind of question ;-) –  rasher Jul 11 at 22:07

1 Answer 1

up vote 8 down vote accepted
+100

Let me start by saying that this problem is probably best solved with a procedural backtracking algorithm, like the one given here. This makes Mathematica a poor choice for tackling it. In fact, judging from this sci.math topic from 1994 the people who originally derived the complete list of PPDIs did it in C.

But since we're here to talk about Mathematica, let's do that. A big efficiency hurdle in your code is the uncompiled use of IntegerDigits in the f function. If you compile that into test2 you should already see some improvements.

However, this doesn't work when the PPDIs are no longer machine-sized integers (around 20 digits on my system). When this happens IntegerDigits can no longer be used in compiled form, so we need another fast way of checking whether a given count of digits corresponds to a PPDI or not.

Here's my attempt at an improved version of your test2:

ClearAll[ValidVectorQ];
With[{
   range = Range[0, 9],
   powerdigits = Reverse /@ IntegerDigits[Range[0, 9]^n, 10, n],
   min = N[10^(n - 1)],
   max = N[10^n],
   powers = Range[0, 9]^N@n
   },
  ValidVectorQ = Compile[
     {{vector, _Integer, 1}},
     If[min <= vector.powers < max,
      Catch[
       Boole@SameQ[
         vector,
         Last /@ Tally[range~Join~(Mod[
                Rest@FoldList[Quotient[#1, 10] + #2 &, 0,
                 If[vector[[Mod[Last[#], 10] + 1]] === 0,
                    Throw[0]; #,
                    #
                    ] &[vector.powerdigits]
                 ],
                10])] - 1
         ]
       ],
      0
      ],
     CompilationTarget -> "C", RuntimeOptions -> "Speed"
     ];
  ];

Instead of using your way to represent PPDIs by the count of their digits, I'm using integer vectors of length 10 whose the $i$th entry indicate how many digits of $i-1$ there are in the PPDI. Thus for 92727 the vector would be {0, 0, 2, 0, 0, 0, 0, 2, 0, 1}.

ValidVectorQ first uses a version of your test1 to check if the PPDI is in the valid size range. It then proceeds to compute all digits of the PPDI recursively with Mod[Rest@FoldList[Quotient[#1, 10] + #2 &, 0, vector.powerdigits], 10]. Once the digits are computed, it computes the corresponding vector with Last /@ Tally[range~Join~digits] - 1 (which I modified from this answer) and checks for equality with the input vector. Note that it doesn't return a boolean but rather 0 or 1 -- this is faster when used in combination with Pick (as explained here).

The rest of my code more or less follows your strategy. It computes the PPDIs of length 21 in roughly 1 minute, and those of length 39 in 80 minutes. For sake of completeness, here it is (it's a bit long, sorry):

n = 39;

PartitionInRangeQ[partition_] :=
  And[
   10^(n - 1) <= partition.Range[10 - Length@partition, 9]^n,
   10^n > Reverse[partition].Range[0, Length@partition - 1]^n
   ];

partitions = 
  Select[Sort /@ IntegerPartitions[n, 10], PartitionInRangeQ];

partitionsByLength = 
  Table[Select[partitions, Length[#] === i &], {i, 1, 10}];

ClearAll[PartitionInRangeQ2];
With[
    {
     min = N[10^(n - 1)],
     max = N[10^n],
     length = #,
     bigpowers = Range[10 - #, 9]^N@n,
     smallpowers = Range[0, # - 1]^N@n
     },
    PartitionInRangeQ2[length] =
     Compile[
      {{partition, _Integer, 1}},
      Boole@And[
        min <= partition.bigpowers,
        partition.smallpowers < max
        ],
      CompilationTarget -> "C", RuntimeOptions -> "Speed"
      ]
    ] & /@ Range[9];

ClearAll[ProcessPartition];
Table[
  With[{length = i},
   ProcessPartition[length] = 
     SelectValidVectors[
       Join @@ (InsertZeros[
           length] /@ (Pick[#, PartitionInRangeQ2[length] /@ #, 1] &[
            Permutations@#]))
       ] &;

   ],
  {i, 1, 9}
  ];
ProcessPartition[10] = SelectValidVectors@Permutations@# &;

ClearAll[InsertZeros];
Table[
  With[{
    insertPositions = 
     Thread@List[1 + # - Range[10 - i]] & /@ 
      Subsets[Range[10], {10 - i}],
    length = i
    },
   InsertZeros[length] = Compile[
      {{partition, _Integer, 1}},
      Insert[partition, 0, #] & /@ insertPositions,
      CompilationTarget -> "C", RuntimeOptions -> "Speed"
      ];
   ],
  {i, 1, 9}
  ];

ClearAll[SelectValidVectors];
SelectValidVectors = Pick[#, ValidVectorQ /@ #, 1] &;

validVectors = 
  Union @@ Flatten[
    Table[ProcessPartition[i] /@ partitionsByLength[[i]], {i, 1, 10}],
     1]

ppdis = validVectors.Range[0, 9]^n
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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