PWC #174

Challenge 1 (Disarium Numbers)

For challenge 1, we are tasked to find the first 19 disarium numbers. These are base-ten positive integers "abc..." such that a^1+b^2+c^3 ... = abc... (I don't know how to type Math in this blog platform, so I use something like Excel formula syntax, with which readers should be familiar)

For example, 518 is a disarium number, because 5+1^2+8^3 = 518. 

There are only 20 disarium numbers, with the 20th being a 20-digit monster. The first 10 are mechanically the single-digit numbers from 0 through 9. The 19th is a 7-digit number around 2.6 million.

I used a brute-force approach of just going through all the numbers from 0 to 3_000_000, and selecting the ones which meet the disarium condition. But a more efficient strategy exists, which I sketch below, but did not implement in most of my scripts. The brute force approach is fast enough for the first 19 items.

More efficient algorithm (but I use brute force for my main submission)

Consider if we want to find all the 2-digit disarium numbers. Rather than cycle through all the 2-digit numbers, we actually just need to cycle through the one-digit numbers 0 through 9. 

Say if a candidate 2-digit number is some 'ab' where a+b^2=10a+b. Then re-writing as an equation in a, we get a = (b^2-b)/9. 

Since 'a' is a unique solution to that equation, and 'a' must be a positive integer, we just need to loop through 1-digit numbers 'b', and find those where 'b^2-b' is divisible by 9. Except for the trivial b=0 and b=1, the only candidate is 9, from which we can construct the disarium number 89 by solving a = (81-9)/9.

Extending to 3-digits, say for candidate 'abc', we have:

a+b^2+c^3=100a+10b+c => a = ( (b^2+c^3)-(10b+c) )/99.

Of course '10b+c' is the 2-digit number 'bc'. 

So to find all 3-digit disarium numbers, we loop through all 2-digit numbers 'bc' checking if the expression ( (b^2+c^3) - 'bc') is divisible by 99, and can then find the unique 'a' for each such case.

And so on. It is easy to generalize to any n-digit number up to n=22. (There is a rigorous proof that there are no Disarium numbers with n greater than 22).

Using this algorithm, to find the first 19 Disarium numbers, we just need to loop from 0 through around 650,000 rather than 0 through 2.6+ million, and have a less expensive check for each number.

Implementation (mostly brute force)

I tried PDL for my Perl 5 script. My script constructs a piddle vector containing the numbers from 0 to 3_000_000, and uses a sub called "get_poly" to find the value of a^1+b^2+c^3... for each such number 'abc...'. This sub is written in Perl/PDL, not in the C-like PP preprocessor language. I use a 'where' PDL expression to filter the piddle to the Disarium numbers x where get_poly(x)==x. I used the 'broadcast_define' routine in PDL::Core to set up broadcasting using a Perl 5/PDL sub. 

I expected this to be blazingly fast, but it was not! In fact, it was very slow indeed, taking around 3 minutes 30 seconds. A simple alternative script using plain Perl 5 loops ran in 8 seconds. My Raku version also using Raku 'for' loops ran in around 1 minute 40 seconds. A Julia script that uses the same matrix algebra as my Perl 5 PDL script ran in under 2 seconds. (I also modified my non-PDL Perl 5 script to try the more efficient algorithm given above, it ran in 1.7 seconds.)

The reason for the inefficiency with PDL is that it is not good at vectorizing functions written in Perl, as opposed to PP which gets translated to C. The solution would be to write my get_poly sub using PP instead with the Inline::Pdlpp module. 

I did attempt this (borrowing shamelessly from a C solution submitted by Laurent Rosenfeld here (See also his blog here) Note for example his elegant approach to computing the length of an integer by using the base-10 logarithm). This was fast as expected, around 3 seconds the first time which includes compile time, and around 0.5 seconds thereafter. 

I also tried a PDL::PP version of the faster algorithm above (not uploaded to github, see appendix below). This ran in 0.32 seconds post-compilation and is my fastest solution in Perl 5. I did also try to implement the faster algorithm in Raku (not uploaded to github), but it was even slower than my brute-force Raku solution, running in around 1 minute 50 seconds. Not sure what's going on there. Update: I was able to improve this, the slowness  was largely due to using lazy sequences x...y instead of straight ranges x..y for for loops. I have provided my Raku code below as an appendix. It runs in around 14 seconds.

Here is my Perl 5 PDL script.

Here is my faster Perl 5 PDL script using PP

Here is my alternative Perl 5 script using vanilla Perl 5 loops.

Here is my Raku script using vanilla loops.

Here is my Julia script using PDL-like matrix algebra.

Here is my Perl 5 non-PDL script using the faster algorithm  

Fast alternative algorithm from another contributor

I noticed this fast algorithm used in the Raku script from Flavio Poletti alias Polettix. It runs in under a second with Raku. It is much more efficient than my fast algorithm above. The logic is to set lower and upper bounds for candidates of different lengths (2-digit, 3-digit, 4-digit, etc.) and recursively bring these closer together, iterating through candidates only when the gap is less than 10. As he suggests, the algorithm could be used to find the big 20th number in reasonable time if translated to a fast compiled implementation with parallelization. I tried running his script to find the 20th number, and gave up while it was scanning through the 17-digit numbers. It slows down for the bigger lengths, not just because of the increase in scale, but because Raku switches to slower bigint arithmetic.

Challenge 2 (rank permutations)

Challenge 2 asks us to create 2 subroutines: (1) "rank2permutation" which given a list 'a' and a positive integer 'n', finds the particular permutation of the items in 'a' that have rank 'n' among all permutations in lexicographic (dictionary) order; and (2) "permutation2rank" which given a list 'a' of items in some order, finds the particular rank that 'a' would have among all permutations of those items in lexicographic order.

There are libraries that make this easy in Perl 5 and Julia, while for Raku, permutations are built-in. For Perl 5 I used the List::Permutor library. For Julia I used the Combinatorics library. (Except for Raku, it's not clear that these libraries particularly use lexicographic order by default. But they all work for the example given in the challenge (for which lexicographic and numeric order are the same)).

Here is my Perl 5 script.

Here is my Raku script.

Here is my Julia script.

Other contributors

Philippe Bricout alias brxfork has submitted a regex-based solution in Perl 5. Amazing again what one can do with regex. 

 

Appendix: my fast PDL+PP version of ch-1

Here is the code for my faster-algorithm version of the challenge 1 solution using PDL and PDL::PP (not uploaded to github: I finished it too late). The inline PP code is ugly, but it works and runs fastest compared to my other scripts for this challenge: 0.32 seconds.

 #-- Listing of alt2-ch-1-pp.pl

#!/usr/bin/env perl

 

#real   0m0.319s

#user   0m0.290s

#sys    0m0.027s

 

#-- PDL solution using inline PP

 

use PDL;

use PDL::NiceSlice;

use PDL::AutoLoader;

use Inline 'Pdlpp';

 

#--find the sequence by checking ints up to 650_000

 

my $pdl = sequence(650_000);

print sequence(10),"\n";

print $pdl->where($pdl->disaurium)->disaurium,"\n";

 

#[0 1 2 3 4 5 6 7 8 9]

#[89 135 175 518 598 1306 1676 2427 2646798]

 

__END__

 

__Pdlpp__

 

#line 1 "my-inline-work"

 

pp_addhdr('

#include "math.h"

');

 

pp_def('disaurium',

    Pars => 'a(); [o]b()',

    Code => pp_line_numbers(__LINE__,q{

        long mya_rw = $a();

        long mya_r=$a();

        long mya_length = ( (mya_r <= 9) ? 1 : (floor(log10($a()))+1) );

        long len_r = mya_length;

        long myb=0;

        long myx=(pow(10,len_r)-1);

        while (mya_rw > 0) {

            myb += pow(mya_rw % 10, mya_length+1);

            mya_rw /= 10;

            mya_length--;

        }

        if ( (((myb-mya_r) % myx)==0) &&

            ((myb - mya_r) > 0) ) {

                $b()= ((myb-mya_r)/myx)*pow(10,len_r)+mya_r;

            }

    }),   

);

Update: Appendix 2: my faster-algorithm Raku script for ch-1

Here is an updated Raku script that uses the faster algorithm I mentioned above. It runs in around 14 seconds on my hardware. Of course, the much faster algorithm from Flavio Poletti runs in less than a second.

 

#!/usr/bin/env raku

#real   0m13.768s
#user   0m13.700s
#sys    0m0.058s


(0 .. 9).map({.say});

(^Inf).grep({&disarium($_)}).head(9).map({say &disarium($_)});

sub disarium ($n) {
    #-- inputs
    my @nstr=$n.comb;
    
    my $retval=0;
    
    #-- calculate a^2+b^3+c^4 ... for $n=abc...
    for (0 .. (@nstr.elems - 1)) -> $i {
        $retval += (@nstr[$i] ** ($i+2));
      }
 
      #check if $retval is divisible by 10^length($n)-1
      my $x = (10 ** @nstr) - 1;
     
     
      #if divisible construct and return disarium number
      if ( ( ($retval - $n) %% $x ) && (($retval-$n) > 0) ) {
          return ( ($retval-$n) / $x ) ~ $n;
      }
      return Nil;
}