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

Way more work than needed - cut your search space to all the non-decreasing sequences of

Ndigits, check these for matching sums after raising to power (pre-compute the latter), if sum sorted non-decreasing matches sequence, you have a hit... – ciao – 2014-07-09T07:05:17.887@rasher I'm afraid that's not sufficiently efficient. For

Ndigits you'll need to sift through`Binomial[N+9,9]`

non-decreasing sequences, which gives you 1.677.106.640 combinations forN= 39. – Teake Nutma – 2014-07-11T08:55:59.527@TeakeNutma: Add parity check, some other rudimentary sieving, you're at a few hundred million checks for

allthe 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 ;-) – ciao – 2014-07-11T22:07:57.293