# Perl Weekly Challenge 147.

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

``````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  #
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  }
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.

``````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

``````EmptyEmptyEmptyEmpty
``````

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  #
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  #
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  #
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