PWC #174

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


Comments

Post a Comment

Popular posts from this blog

PWC 258

PWC 258 Last week (#257) I finished challenge 2 and my blog post for challenge 257 late last week, so had only submitted challenge 1. Here is my blog post for PWC 257 , including the full code for my challenge 2 solution. This week, I have done both challenges in perl 4.019 on DOSBOX. Both scripts also work fine on the latest Perl 5.38. I also did the challenges in Raku, where they are one-liners.  Challenge 1 (Count Even Digits Number) We are given a list of integers, and asked to count the number of elements whose base-10 representation has an even number of digits. In Raku , this is a one-liner. Here it is: raku -e 'say +@*ARGS.grep({$_.chars %% 2})' $@ It could be a one-liner in Perl 4 too, but since I run it on DOSBOX, my script would not be a bash shell script. Here is the short subroutine that does the job.  sub count_even_digits {  scalar( grep( ((length($_) % 2) == 0), @_) ); } The full code is here . Challenge 2 (Sum of values) We are given an integer k  and a list in
Read more

PWC 249

PWC 249 Challenge 1 (Equal pairs) We are given a list of integers with an even number of elements. We are asked to check if the list can be completely divided into sub-lists of size 2 with each sub-list containing the same integer repeated, like say (2,2). If so, we are asked to return the list of such sub-lists. This is straightforward in Perl 5. I use modern conveniences this time, though it could be done using only primitive Perl 4 syntax. Here is my subroutine. It returns a list of array references. local *equal_pairs = sub {     #-- return () if input is of odd length     ((scalar(@_) % 2) > 0) && (return ());          #-- count number of occurrences of each element using a hash     local %_ = ();     map {$_{$_}++} @_;          #-- if any element occurs an odd number of times, return ()     (grep {($_{$_} % 2) > 0} keys %_) && (return ());          #-- loop thru %_ keys, returning pairs of such keys     map {         my $k=$_;         my $num_pairs=( ($_{$k}
Read more

PWC 253

PWC 253 I do both challenges this week using old software on DOSBOX: perl 4.019, and for comparison, also python 1.4 beta. My solutions also work on the current versions of Perl 5 (5.38) and Python 2 (2.7.18) respectively. Challenge 1 (Split Strings) Links to code: perl4.019 python1.4b We are given a list of strings, and a single character say sep . We are asked to return the list, with the strings split into component substrings using sep as a separator. We should not include any empty strings in the output. In perl 4, I first stick a backslash in front of the separator character, so that it is escaped if entered in a regex. This is because the perl split function takes a regex as its first argument indicating the separator. If we have a character like '$' or ',', we want to escape it so that it is read literally when parsing the regex. Here is the key subroutine:   sub split_strings {  local ($separator, @words)=@_;  local (@retval)=(); #array to be returned  #  #--
Read more