# Perl Weekly Challenge 174.

``````Submitted by: Mohammad S Anwar
Write a script to generate first 19 Disarium Numbers.

A disarium number is an integer where the sum of each digit
raised to the power of its position in the number, is equal to
the number.

For example,

518 is a disarium number as (5 ** 1) + (1 ** 2) + (8 ** 3) =>
5 + 1 + 512 => 518
``````

I try a straightforward solution examining numbers one by one and testing them until I find enough Disarium numbers. I use PDL for conciseness.

``````perl  -MPDL -E '
\$N=shift; \$n=\$i=0; while(\$n<\$N){\$e=(\$p=long [split "", \$i])->sequence+1;
++\$n, say "\$n: \$i" if (\$p**\$e)->sumover==\$i; ++\$i }' 19
``````

Results:

``````1: 0
2: 1
3: 2
4: 3
5: 4
6: 5
7: 6
8: 7
9: 8
10: 9
11: 89
12: 135
13: 175
14: 518
15: 598
16: 1306
17: 1676
18: 2427
19: 2646798
``````

It seems to work fine and is fast for the first 18 desarium numbers, but the 19-th took about 2.5 minutes in my laptop. I’m not sure I would have the patience for waiting for the 20-th. The full code is

`````` 1  # Perl weekly challenge 174
2  # Task 1: Disarium numbers. Brute force
3  #
5  use v5.12;
6  use PDL;
7  say("Usage: \$0 N\nto obtain the first N Disarium numbers"), exit unless @ARGV;
8  my \$how_many=shift;
9  die "Expected a positive number" unless \$how_many>0;
10  my (\$counter, \$test)=(0, 0);
11  while(\$counter<\$how_many){
12      my \$exponents=(my \$digits=long [split "", \$test])->sequence+1;
13      ++\$counter, say "\$counter: \$test"
14          if (\$digits**\$exponents)->sumover==\$test;
15      ++\$test
16  }
``````

Example:

``````./ch-1.pl 19
``````

Results:

``````1: 0
2: 1
3: 2
4: 3
5: 4
6: 5
7: 6
8: 7
9: 8
10: 9
11: 89
12: 135
13: 175
14: 518
15: 598
16: 1306
17: 1676
18: 2427
19: 2646798
``````

After a twitter exchange with Flavio Poletti `@polettix` I learned that one can do better than a brute force approach. I saw his Raku solution PWC174 - Disarium Numbers - ETOOBUSY not quite understanding it, but it inspired the following solution.

First, note that a number given by a sequence of N digits, y1y2…yN with a positive first digit is Disarium if p(y1y2…yN)=d(y1y2…yN), where p(y1y2…yN)=∑n=1N ynn is the sum of powers of its digits and d(y1y2…yN)=∑ 10N-n yn is its decimal value.

Notice that p(x1…xMyM+1…yN)=p(x1…xM)+p(01…0MyM+1…yN) and d(x1…xMyM+1…yN)=10N-M d(x1…xM)+d(yM+1…yN), where 01…0M represents M zeroes, one after another. Thus, d(x1…xM)=(p(01…0MyM+1…yN)-d(yM+1…yN)+p(x1…xM))/10N-M.

The smallest possible value for p(x1…xM) is 1, corresponding to 1102…0M, while the largest possible value is ∑n=1M 9n=(9M+1-9)/8. Thus, (p(01…0MyM+1…yN)-d(yM+1…yN)+1)/10N-M ≤ d(x1…xM) ≤ (p(01…0MyM+1…yN)-d(yM+1…yN)+(9M+1-9)/8)/10N-M.

The strategy is thus, choose a value for yN, and find the bounds for the rest of the digits x1…xN-1. If the bounds are sufficiently tight, then explore all intermediate numbers looking for Disarium numbers. If the bounds are too far apart, then choose also a value for yN-1 and find bounds for x1…xN-2. Recourse until the bounds are tight enough and then explore exhaustively the given range. For each M explored, starting downwards from N-1 and stopping when tight enough bounds are found, Choose all possible digits 0…9 for YM+1.

These ideas lead to the following recursive code.

`````` 1  # Perl weekly challenge 174
2  # Task 1: Disarium numbers. Not so brute force
3  #
5  use v5.36;
6  use List::Util qw(sum min);
7  use Memoize;
8  memoize "power";
9  die "Usage: \$0 N\nto obtain the first N Disarium numbers\n" unless @ARGV;
10  my \$want=shift;
11  die "There are not that many Disarium numbers\n" if \$want > 20;
12  die "We are not that patient today\n" if \$want==20;
13  my \$count=1;
14  say \$count++, ": \$_" for(0..min 9, \$want);
15  for my \$N(2..20){ # for each possible length
16      search(\$N, \$N, "", 0, 0);
17  }
18  sub search (\$N, \$M, \$rightmost, \$p_old, \$d_old){
19      --\$M;
20      my \$minmin=power(10,\$M-1);
21      for my \$right(0..9){ # for each possible rightmost digit
22          last if \$count > \$want;
23          my \$p_new=power(\$right, \$M+1)+\$p_old;
24          my \$d_new=\$right*power(10, \$N-\$M-1)+\$d_old;
25          my \$min=(\$p_new-\$d_new+1)/power(10, \$N-\$M);
26          my \$max=(\$p_new-\$d_new+(power(9,\$M+1)-9)/8)/power(10, \$N-\$M);
27          \$min=\$minmin if \$min<\$minmin;
28          next if \$max<=\$min;
29          my \$newright="\$right\$rightmost";
30  	if(\$max-\$min < 10){
31  	    for(\$min..\$max){
32                  my \$candidate="\$_\$newright";
33                  say \$count++, ": \$candidate" if disarium(\$candidate);
34              }
35          } else {
36              search(\$N, \$M, \$newright, \$p_new, \$d_new);
37          }
38      }
39  }
40  sub power(\$x, \$y){ \$x**\$y }
41  sub disarium(\$candidate){
42      my @digits=split "", \$candidate;
43      my \$powersum=sum map {\$digits[\$_-1]**\$_} (1..@digits);
44      return \$powersum==\$candidate;
45  }
``````

I `memoize` the `sub power` to avoid recalculating the same powers again and again. Line 10 deals with the 1 digit case. The recursive subroutine `search` takes as parameters the number of digits N, the numbers of unknown digits M, a string with the rightmost N-M candidate digits yM+1…YN, and the corresponding values of p(yM+1…YN) and d(yM+1…YN). Then it decreases M and starts exploring values of the new yM+1, calling itself recursively if necessary.

Example:

``````./ch-1a.pl 19
``````

Results:

``````1: 0
2: 1
3: 2
4: 3
5: 4
6: 5
7: 6
8: 7
9: 8
10: 9
11: 89
12: 135
13: 175
14: 518
15: 598
16: 1306
17: 1676
18: 2427
19: 2646798
``````

So it works. The above run took 0.026s instead of 2.5m for the previous code, so it is about six thousand times faster.

To get the last (20-th) Disarium number I would require `bigint` and lots of patience. Using `bigint` and `Time::HiRes` I timed the time required for a few lengths (up to 14) and extrapolated. My estimate is that it would take about 142 hours to get the 20-th Disarium number.

Curiously, the brute force solution but without using PDL is only about ten times slower for the first 19 Disarium numbers than the last, somewhat involved program, so it seems that in this case PDL is not so helpful.

``````Submitted by: Mohammad S Anwar
You are given a list of integers with no duplicates, e.g. [0,
1, 2].

Write two functions, permutation2rank() which will take the
list and determine its rank (starting at 0) in the set of
possible permutations arranged in lexicographic order, and
rank2permutation() which will take the list and a rank number
and produce just that permutation.

algorithm.

Given the list [0, 1, 2] the ordered permutations are:

0: [0, 1, 2]
1: [0, 2, 1]
2: [1, 0, 2]
3: [1, 2, 0]
4: [2, 0, 1]
5: [2, 1, 0]
and therefore:

permutation2rank([1, 0, 2]) = 2

rank2permutation([0, 1, 2], 1) = [0, 2, 1]
``````

A straightforward solution may be easily constructed by producing all permutations in order, which yields the following 3-liner:

``````perl -MAlgorithm::Combinatorics=permutations -MPDL -MPDL::NiceSlice -E '
\$p=pdl "[".shift."]"; @o=\$p->qsort->list; \$a=pdl[permutations(\@o)];
for(0..\$a->dim(1)-1){say("rank=\$_"), last if all(\$a(,\$_)==\$p);}
say "\$_->",\$a(:,(\$_)) foreach(@ARGV);
' "1 0 2" 0  2  4 1 3 5
``````

Results:

``````rank=2
0->[0 1 2]
2->[1 0 2]
4->[2 0 1]
1->[0 2 1]
3->[1 2 0]
5->[2 1 0]
``````

Similar example with a different set of numbers:

``````perl -MAlgorithm::Combinatorics=permutations -MPDL -MPDL::NiceSlice -E '
\$p=pdl "[".shift."]"; @o=\$p->qsort->list; \$a=pdl[permutations(\@o)];
for(0..\$a->dim(1)-1){say("rank=\$_"), last if all(\$a(,\$_)==\$p);}
say "\$_->",\$a(:,(\$_)) foreach(@ARGV);
' "10 5 15" 0  2  4 1 3 5
``````

Results:

``````rank=2
0->[5 10 15]
2->[10 5 15]
4->[15 5 10]
1->[5 15 10]
3->[10 15 5]
5->[15 10 5]
``````

It works but for larger sets it becomes extremely inefficient, as the number of permutations grows as the factorial of the number of elements. Thus, for the full code I change to an algorithm based on the page The rank of a permutation, TryAlgo.

Consider a permutation P0P1…Pn-1 of the ordered set O0O1…On-1. There are (n-1)! permutations that start with O0, (n-2)! that start with O1 and so on. Thus, the rank r(P0P1…Pn-1) is simply r0×(n-1)!+r(P1…Pn-1), where r0 is the sequence number of P0 in O0O1…On-1, counting up from 0, and r(P1…Pn-1) is the rank of the permutation P1…Pn-1 of the ordered set O0O1…On-1 from which we have removed P0. We can proceed recursively to obtain finally r(P0P1…Pn-1)=r0×(n-1)!+r1×(n-2)!+…rn-1×0!, where ri is the position of Pi in the set O0O1…On-1 from which we have removed P0, P1…Pi-1. Notice that rn-1=0 as there is no choice left for the last element.

Conversely, given a rank r, we can obtain r0 by dividing it by (n-1)!, we can obtain r1 by dividing the residue by (n-2)!, and so on. Then, we choose P0 as the r0-th element of the set O0O1…On-1, we remove P0 from the set and choose P1 as the r1-th element of the remaining set, and so on. This algorithm is embodied in the following code:

`````` 1  # Perl weekly challenge 174
2  # Task 2: Permutation ranking
3  #
5  use v5.36;
6  use PDL;
7  use PDL::NiceSlice;
8  my \$usage=<<"END";
9  Usage: \$0 "P1 P2...Pn" R1 R2...
10  to obtain rank of permutation P1 P2...Pn
11  and the rank-Ri permutations. Note quotations.
12  END
13  say(\$usage), exit unless @ARGV>0;
14  my \$permutation=long "[".shift."]";
15  my \$size=\$permutation->nelem;
16  die "Elements should be distinct" unless \$permutation->uniq->nelem==\$size;
17  my @ordered=\$permutation->qsort->list;
18  my %element2index = map { (\$ordered[\$_], \$_) } (0..\$size-1);
19  my \$permuted_indices=long [@element2index{\$permutation->list}];
20  my \$factorials=factorials(\$size);
21  say "permutation2rank(\$permutation)=", permutation2rank(\$permuted_indices);
22  say "rank2permutation(\$permutation, \$_)=", rank2permutation(\$_) for(@ARGV);
23
24  sub permutation2rank(\$permutation){ # ranks a permutation of 0...N-1
25      return \$factorials->inner(ranks(\$permutation));
26  }
27  sub ranks(\$permutation){
28      my \$cmp=\$permutation(*1)>\$permutation; #P_i>P_j
29      \$cmp->inner(\$cmp->xvals>=\$cmp->yvals); #r_i=sum_{j>=i}(P_i>P_j)
30  }
31  sub factorials(\$size){
32       my \$f=1;
33       (long [1, map {\$f=\$_*\$f} (1..\$size-1)])->(-1:0);
34  }
35  sub rank2permutation(\$rank){
36      my @indices=map {my \$index=floor(\$rank/\$_); \$rank%=\$_; \$index} \$factorials->list;
37      my @copy=@ordered;
38      long [ map {splice @copy, \$_, 1} @indices];
39  }
40
``````

The input to the program is a list of distinct numbers between quotation marks, that define the permutation followed by a list of ranks to examine. I use PDL to simplify some expressions. In line 17 I produce the ordered set O0O1…O{n-1}, which in lines 18 and 19 I map to the simpler set 0,1…n-1, for which the indices and correponding values coincide. Line 20 prepares a vector of decreasing factorials by calling the subroutine `factorials`. The subroutine `ranks` calculates a vector of partial ranks ri by counting how many Pj’s obey Pj>Pi, with i≥j. To that end, in line 28 we use a dummy index to build a matrix with entries Pi>Pj at row i column j. Then, we use an inner product in line 29 to select and sum its upper right triangle, yielding ri=∑j>=i(Pi>Pj). The subroutine `rank2permutation` maps a rank r to a permutation by iteratively obtaining an array of indices ri=r/fi updating the rank r→r%fi. I use the `splice` routine in line 38 to chose and remove the ri-th term Pi of the set, and return a vector with P0P1…Pn-1.

Example:

``````./ch-2.pl "1 0 2" 0 1 2 3 4 5
``````

Results:

``````permutation2rank([1 0 2])=2
rank2permutation([1 0 2], 0)=[0 1 2]
rank2permutation([1 0 2], 1)=[0 2 1]
rank2permutation([1 0 2], 2)=[1 0 2]
rank2permutation([1 0 2], 3)=[1 2 0]
rank2permutation([1 0 2], 4)=[2 0 1]
rank2permutation([1 0 2], 5)=[2 1 0]
``````

Another example, similar to above but with different numbers:

``````./ch-2.pl "10 9 20" 0 1 2 3 4 5
``````

Results:

``````permutation2rank([10 9 20])=2
rank2permutation([10 9 20], 0)=[9 10 20]
rank2permutation([10 9 20], 1)=[9 20 10]
rank2permutation([10 9 20], 2)=[10 9 20]
rank2permutation([10 9 20], 3)=[10 20 9]
rank2permutation([10 9 20], 4)=[20 9 10]
rank2permutation([10 9 20], 5)=[20 10 9]
``````

The run-time of this algorithm is proportional to the square of the number of elements, so it is much faster than the oneliner above.

Notice that by replacing line 17 by

``````my @ordered=sort {\$a cmp \$b} \$permutation->list;
``````

I would obtain a lexicographic order instead of a numerical order (see task 2a).

Example:

``````echo "Numerical order"
./ch-2.pl "22 111 3" 0 1 2 3 4 5
echo; echo "Lexicographical order (with modified line 17)"
./ch-2a.pl "22 111 3" 0 1 2 3 4 5
``````

Results:

``````Numerical order
permutation2rank([22 111 3])=3
rank2permutation([22 111 3], 0)=[3 22 111]
rank2permutation([22 111 3], 1)=[3 111 22]
rank2permutation([22 111 3], 2)=[22 3 111]
rank2permutation([22 111 3], 3)=[22 111 3]
rank2permutation([22 111 3], 4)=[111 3 22]
rank2permutation([22 111 3], 5)=[111 22 3]

Lexicographical order (with modified line 17)
permutation2rank([22 111 3])=2
rank2permutation([22 111 3], 0)=[111 22 3]
rank2permutation([22 111 3], 1)=[111 3 22]
rank2permutation([22 111 3], 2)=[22 111 3]
rank2permutation([22 111 3], 3)=[22 3 111]
rank2permutation([22 111 3], 4)=[3 111 22]
rank2permutation([22 111 3], 5)=[3 22 111]
``````

Furthermore, I can also generalize to non-numerical inputs by replacing some ndarrays by perl arrays (see task 2b).

``````./ch-2b.pl "BB AAA C" 0 1 2 3 4 5
``````

Results:

``````permutation2rank([BB AAA C])=2
rank2permutation([BB AAA C], 0)=[AAA BB C]
rank2permutation([BB AAA C], 1)=[AAA C BB]
rank2permutation([BB AAA C], 2)=[BB AAA C]
rank2permutation([BB AAA C], 3)=[BB C AAA]
rank2permutation([BB AAA C], 4)=[C AAA BB]
rank2permutation([BB AAA C], 5)=[C BB AAA]
``````
Written on July 18, 2022