Perl Weekly Challenge 147.

My solutions (task 1 and task 2 ) to the The Weekly Challenge - 147.

Task 1: Truncatable Prime

Submitted by: Mohammad S Anwar
Write a script to generate first 20 left-truncatable prime
numbers in base 10.

In number theory, a left-truncatable prime is a prime number
which, in a given base, contains no 0, and if the leading left
digit is successively removed, then all resulting numbers are
primes.

Example
9137 is one such left-truncatable prime since 9137, 137, 37
and 7 are all prime numbers.

One could solve this problem by doing a kind of Eratosthenes sieve, first eliminating non-primes, and then eliminating non-truncable primes. This may be done with PDL in a three liner. I’ll leave the explanation for the comments of the full code below.

perl -d -MPDL -MPDL::NiceSlice -E ' $N=200; $s=ones($N); $s(0:1).=0;
$s($_**2:-1:$_).=0 foreach(2..sqrt($N-1)); for($n=10; $n<$N; $n*=10){
$s->reshape($n,$N/$n); $s&=$s(:,0); $s(:,10:-1:10).=0 if $s->dim(1)>10;}
$s->reshape($N); say $p=$s->xvals->where($s)->(0:19);'

Results:

[2 3 5 7 13 17 23 37 43 47 53 67 73 83 97 113 137 167 173 197]

The full code is the following:

 1  # Perl weekly challenge 147
 2  # Task 1: Truncatable prime
 3  #
 4  # See https://wlmb.github.io/2022/01/10/PWC147/#task-1-truncatable-prime
 5  use v5.12;
 6  use warnings;
 7  use PDL;
 8  use PDL::NiceSlice;
 9  use POSIX (); # don't import to avoid name collisions with PDL
10  use Text::Wrap qw(wrap $columns $break);
11
12  die "Usage: ./ch-1.pl size_of_sieve number_of_truncatable_primes [base]\n"
13      unless @ARGV>=2;
14  my ($size, $wanted, $base)=@ARGV;
15  $base//=10; # decimal numbers by default
16  $size=$base**POSIX::ceil(log($size)/log($base)); # Force power of base;
17  my $sieve=ones($size); #
18  $sieve(0:1).=0; # 0 and 1 are not primes
19  # find primes with Eratosthenes sieve
20  $sieve($_**2:-1:$_).=0 foreach(2..sqrt($size-1));
21  # Remove non-truncatable primes
22  for(my $n=$base; $n<$size; $n*=$base){
23      # Reshape sieve as rectangle. The first row all log_base(n) digits
24      $sieve->reshape($n,$size/$n);
25      # Remove from the remaining rows those numbers which don't
26      # correspond to a truncatable prime in the first row
27      $sieve &= $sieve(:,0);
28      # From every tenth row remove numbers that would begin in 0 if truncated
29      $sieve(:,10:-1:10).=0 if $sieve->dim(1)>10;
30  }
31  $sieve->reshape($size); # return to a 1D sieve
32  # The desired primes correspond to the surviving ones in the sieve
33  my $truncatable_primes=$sieve->xvals->where($sieve);
34  my $found=$truncatable_primes->nelem; # truncatable primes actually found
35  say("Didn't find enough ($wanted) primes, please increase size of sieve"),
36      $wanted=$found
37      unless $found >= $wanted;
38  $columns=62; $break=qr/\s/;
39  say wrap("", "    ", "The first $wanted truncatable primes are: ",
40      join ", ", $truncatable_primes(0:$wanted-1)->list);

The idea is to first run the sieve of Erastothenes producing a list of 1’s corresponding to primes and 0’s corresponding to non-primes, as

0  1  2  3  4  5  6  7  8  9  10  11  12  13  14 ...
-----------------------------------------------------
0  0  1  1  0  1  0  1  0  0   0   1   0   1   0 ...

I fold the list as a table of 10 columns, as

    | 0  1  2  3  4  5  6  7  8  9
----------------------------------
 00 | 0  0  1  1  0  1  0  1  0  0
 10 | 0  1  0  1  0  0  0  1  0  1
 20 | 0  0  0  1  0  0  0  0  0  1

An and with the first row removes primes that don’t yield primes when truncated to one digit,

    | 0  1  2  3  4  5  6  7  8  9
----------------------------------
 00 | 0  0  1  1  0  1  0  1  0  0
 10 | 0  0  0  1  0  0  0  1  0  0
 20 | 0  0  0  1  0  0  0  0  0  0

The row number 10 (the eleventh row) would correspond to numbers such as 100, 101, 102… which would have a leading 0 if truncated to two digits and thus are non-truncatable and should be removed. The same is true for every tenth row afterwards.

The procedure above is repeated for a rectangular array with 100 columns

    | 0  1  2  ...  10  11  12  ...  20  21  ...  99
----------------------------------------------------
000 | ...
100 |
200 |

and then 1000, etc.

Example:

./ch-1.pl 100 20

Results:

Didn't find enough (20) primes, please increase size of sieve
The first 15 truncatable primes are: 2, 3, 5, 7, 13, 17, 23,
    37, 43, 47, 53, 67, 73, 83, 97

A sieve of size 100 was not large enough, so I try again with a larger one.

./ch-1.pl 1000 20

Results:

The first 20 truncatable primes are: 2, 3, 5, 7, 13, 17, 23,
    37, 43, 47, 53, 67, 73, 83, 97, 113, 137, 167, 173, 197

I can explore truncation in other bases, for example, in base 8:

./ch-1.pl 1000 20 8

Results:

The first 20 truncatable primes are: 2, 3, 5, 7, 11, 13, 19,
    23, 29, 31, 37, 43, 47, 53, 59, 61, 67, 71, 101, 107

Note that 17 is truncatable in base 10, as 7 is prime, but it is not truncable in base 8, as 17=021 which truncates to 1, while 19 is not truncatable in base 10 as 9 is not prime, but it is truncatable in base 8 as 19=023 which truncates to 3, which is prime.

Task 2: Pentagon Numbers

Submitted by: Mohammad S Anwar
Write a sript to find the first pair of Pentagon Numbers whose
sum and difference are also a Pentagon Number.

Pentagon numbers can be defined as P(n) = n(3n - 1)/2.

Example
The first 10 Pentagon Numbers are:
1, 5, 12, 22, 35, 51, 70, 92, 117 and 145.

P(4) + P(7) = 22 + 70 = 92 = P(8)
but
P(4) - P(7) = |22 - 70| = 48 is not a Pentagon Number.

I made several attempts to solve this task. One didn’t work, one was slow at first (then got faster) and the other two were relatively fast and compact.

First. PDL one-liner

My first attempt was to use PDL to generate pentagon numbers, add them and substract them and check if the agree with some other pentagon number. If I denote the i-th pentagon number as pi, then if pi-pj=pk and pi+pj=pl, then pi=pj+pk and pi+pj=2pj+pk. Thus, I generate a vector p with components pn=n(3n-1)/2 and use them to generate a cube T with components Tjki=pj+pk-pi. I want those indices for which Tjki=0. Similarly, I build another cube U with components Ujkl=2pj+pk-pl and look for the indices for which Ujkl=0. Finally, I can construct the solution from those pairs jk common to both sets. I can summarize both conditions by building a Boolean hypercube H with components Hjkil=(pj+pk-pi==0)&(2pj+pk-pl=0) and look for the indices jkil of the true values. From those indices and the vector p of pentagon numbers I can immediately obtain the solution. This is readily translated to a very compact PDL code.

perl -MPDL -MPDL::NiceSlice -E '
$N=20; $n=zeroes($N)->xvals+1; $p=$n*(3*$n-1)/2;
$x=whichND(($p+$p(*1)-$p(*1,*1)==0)&(2*$p+$p(*1)-$p(*1,*1,*1)==0));
($pj,$pk)=map {$p->index($x(($_)))} (0,1); say "\n", $pj+$pk,$pj, $pk, 2*$pj+$pk'

The interesting line is, I believe, $x=whichND(($p+$p(*1)-$p(*1,*1)==0)&(2*$p+$p(*1)-$p(*1,*1,*1)==0)); First the 3D cube $T=$p+$p(*1)-$p(*1,*1) with components $T($j,$k,$i)=$p($j)+$p($k)-$p($i) is built from a vector $p containing the first pentagonal numbers. To that end I use the capability of PDL of transforming vectors to matrices and data cubes by adding leading dummy indices (the *1’s above), and of adding vectors to matrices and data cubes by implicitly adding trailing dummy indices. Comparing to zero I get a cube of Boolean values. Similarly, I make an hypercube with the expression 2*$p+$p(*1)-$p(*1,*1,*1)==0, as the last term has three dummy indices. I combine both with the & operator. Finally, the whichND function yields an ndarray each of whose entries is a list of indices [$j, ~$k, ~$i, ~$l] of pentagon numbers that obey the desired conditions.

Unfortunately, the code above didn’t work. When I ran it I got

Empty[0]Empty[0]Empty[0]Empty[0]

meaning I found no solution. I tried using larger values of $N, but very soon the memory was filled by the 4D hypercube and I had to kill the process. :(

Second: Brute force

Thus I have to change tack. Below I follow a brute force approach. I generate pairs of pentagon numbers in a double loop and stop when I find one that complies with the conditions. To that end I made a routine to identify pentagon numbers. Completing squares, we notice that if p=n(3n-1)/2, then 24p+1=36n2-12n+1=(6n-1)2. Thus, if p is a pentagon number 24p+1 is a perfect square. Furthermore, its square root is s=6n-1, so that s=5 mod 6, and we may easily obtain n=(s+1)/6.

 1  # Perl weekly challenge 147
 2  # Task 2: pentagon numbers
 3  #
 4  # See https://wlmb.github.io/2022/01/10/PWC147/#task-2-pentagon-numbers
 5  use v5.12;
 6  use warnings;
 7  use POSIX qw(floor);
 8  use Time::HiRes qw(time);
 9
10  die "Usage: ./ch-2.pl largest_index\n" unless @ARGV==1;
11  my $N=shift;
12  my $start=time();
13  use integer;
14  J:
15      foreach my $j(2..$N){
16          my $p=$j*(3*$j-1)/2;
17          foreach my $k(1..$j-1){
18              my $q=$k*(3*$k-1)/2;
19              say("p$j=$p\np$k=$q\np$j-p$k=", $p-$q, "=p", index_of($p-$q),
20                    "\np$j+p$k=", $p+$q, "=p", index_of($p+$q)),
21                  last J if pentagonal($q+$p) && pentagonal($p-$q);
22          }
23  }
24  no integer; # don't truncate time
25  say "Time: ", time()-$start;
26  use integer;
27  sub pentagonal {
28      my $p=24*shift()+1;
29      my $s=floor(sqrt($p));
30      return $s**2==$p && $s%6==5;
31  }
32  sub index_of {
33      my $p=24*shift()+1;
34      my $s=sqrt($p);
35      return ($s+1)/6;
36  }
37

./ch-2.pl 100

Results:

Time: 0.00250792503356934

As expected, 100 is too small and produced no result.

./ch-2.pl 500

Results:

Time: 0.064255952835083

$N=500 is still too small, and the time grows with the second power of the size.

./ch-2.pl 2500

Results:

p2167=7042750
p1020=1560090
p2167-p1020=5482660=p1912
p2167+p1020=8602840=p2395
Time: 1.19663405418396

It worked! The sum and difference of the 2167-th and the 1020-th pentagon numbers yield the 2395-th and the 1912-th pentagon number respectively. (My first attempt was much slower, but I had used bigint instead of integer.) I guess the program would have been faster had I checked pi=pj+pk and pl=2pj+pk instead of checking pk=pi-pj and pl=pi+pj.

Third: Remove inner loop

I guess the program can be made faster if I only do explicitly the outer loop and use PDL for the inner loop.

 1  # Perl weekly challenge 147
 2  # Task 2: pentagon numbers
 3  #
 4  # See https://wlmb.github.io/2022/01/10/PWC147/#task-2-pentagon-numbers
 5  use v5.12;
 6  use warnings;
 7  use Time::HiRes qw(time);
 8  use PDL;
 9  use PDL::NiceSlice;
10
11  die "Usage: ./ch-2a.pl largest_index\n" unless @ARGV==1;
12  my $N=shift;
13  my $start=time();
14  my $n=zeroes(long, $N)->xvals+1;
15  my $p=$n*(3*$n-1)/2;
16  for my $i (1..$p->nelem){
17      my $pi=$p(($i-1));
18      my $pass=which(pentagonal($pi+$p) & pentagonal($pi-$p));
19      next unless $pass->nelem;
20      my $j=$pass((0))+1;
21      my $pj=$p(($j-1));
22      my $s=$pi+$pj;
23      my $d=$pi-$pj;
24      my ($k, $l)=map {index_of($_)} ($d, $s);
25      say "p$i=$pi\np$j=$pj\np$i-p$j=$d=p$k\np$i+p$j=$s=p$l";
26      last;
27  }
28  say "Time: ", time()-$start;
29  sub pentagonal {
30      my $p=shift;
31      my $p241=24*$p+1;
32      my $sp241=$p241->sqrt;
33      return (($p>0)&($sp241**2==$p241) & ($sp241%6==5));
34  }
35  sub index_of {
36      my $p=24*shift()+1;
37      my $s=sqrt($p);
38      return ($s+1)/6;
39  }

I run it as

./ch-2a.pl 2500

Results:

p2167=7042750
p1020=1560090
p2167-p1020=5482660=p1912
p2167+p1020=8602840=p2395
Time: 0.912932872772217

It worked too and seems faster, but only by 30%. I had to fiddle the indices to count the pentagon numbers up from one and not from zero.

Fourth: Remove outer loop

I can also get rid of the outer loop using PDL.

 1  # Perl weekly challenge 147
 2  # Task 2: pentagon numbers
 3  #
 4  # See https://wlmb.github.io/2022/01/10/PWC147/#task-2-pentagon-numbers
 5  use v5.12;
 6  use warnings;
 7  use Time::HiRes qw(time);
 8  use PDL;
 9  use PDL::NiceSlice;
10
11  die "Usage: ./ch-2a.pl largest_index\n" unless @ARGV==1;
12  my $N=shift;
13  my $start=time();
14  my $n=zeroes(long, $N)->xvals+1;
15  my $p=$n*(3*$n-1)/2;
16  my $pass=whichND(pentagonal($p+$p(*1)) & pentagonal($p-$p(*1)));
17  die "No solution found. Try to increase the largest_index" unless $pass->dim(1)>0;
18  my $ij=$pass(:,(0))+1;
19  my ($pi, $pj)=map {$p(($_-1))} (my ($i, $j)=map {$ij(($_))} (0,1));
20  my ($s, $d)=($pi+$pj, $pi-$pj);
21  my ($k, $l)=map {index_of($_)} ($d, $s);
22  say "p$i=$pi\np$j=$pj\np$i-p$j=$d=p$k\np$i+p$j=$s=p$l";
23  say "Time: ", time()-$start;
24  sub pentagonal {
25      my $p=shift;
26      my $p241=24*$p+1;
27      my $sp241=$p241->sqrt;
28      return (($p>0)&($sp241**2==$p241) & ($sp241%6==5));
29  }
30  sub index_of {
31      my $p=24*shift()+1;
32      my $s=sqrt($p);
33      return ($s+1)/6;
34  }

Notice that I used again the trick of adding dummy dimensions (line 16), but I only have 2D ndarrays instead of 4D as in the first attempt.

Example:

./ch-2b.pl 2500

Results:

p2167=7042750
p1020=1560090
p2167-p1020=5482660=p1912
p2167+p1020=8602840=p2395
Time: 0.635129928588867

The speed increased another 40%. Overall, using PDL duplicated the speed. Furthermore, the code is somewhat cleaner (this is a subjective statement, of course). A problem with the PDL solutions is that they require more memory and they are exhaustive, i.e., they don’t stop at the first solution.

Fifth: One-liner

Finally, the last solution can be distilled into a one-liner (actually, three-liner):

perl -MPDL -MPDL::NiceSlice -E '$n=zeroes(long,2500)->xvals+1;
$p=$n*(3*$n-1)/2; ($i,$j)=whichND(pent($p+$p(*1))&pent($p-$p(*1)))->(:,0)->list;
($pi,$pj)=($p(($i)),$p(($j))); say "pi=$pi, pj=$pj, pi-pj=", $pi-$pj, " pi+pj=", $pi+$pj;
sub pent {$S=($P=24*shift()+1)->sqrt; $P>0&($S**2==$P)&($S%6==5)}'

The results are obtained again within a fraction of a second:

pi=7042750, pj=1560090, pi-pj=5482660 pi+pj=8602840

The main difference between this program and the first (failed) attempt is that here we use an algebraic test for pentagon-ess instead of determining membership in a given set of pentagon numbers through comparisons.

Written on January 10, 2022