# 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 *p _{i}*, then if

*p*and

_{i}-p_{j}=p_{k}*p*, then

_{i}+p_{j}=p_{l}*p*and

_{i}=p_{j}+p_{k}*p*. Thus, I generate a vector

_{i}+p_{j}=2p_{j}+p_{k}*p*with components

*p*and use them to generate a cube

_{n}=n(3n-1)/2*T*with components

*T*. I want those indices for which

_{jki}=p_{j}+p_{k}-p_{i}*T*. Similarly, I build another cube

_{jki}=0*U*with components

*U*and look for the indices for which

_{jkl}=2p_{j}+p_{k}-p_{l}*U*. Finally, I can construct the solution from those pairs

_{jkl}=0*jk*common to both sets. I can summarize both conditions by building a Boolean hypercube

*H*with components

*H*and look for the indices

_{jkil}=(p_{j}+p_{k}-p_{i}==0)&(2p_{j}+p_{k}-p_{l}=0)*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=36n ^{2}-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
*p _{i}=p_{j}+p_{k}* and

*p*instead of checking

_{l}=2p_{j}+p_{k}*p*and

_{k}=p_{i}-p_{j}*p*.

_{l}=p_{i}+p_{j}## 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.