Perl Weekly Challenge 167.
My solutions (task 1 and task 2 ) to the The Weekly Challenge - 167.
Task 1: Circular Prime
Submitted by: Mohammad S Anwar
Write a script to find out first 10 circular primes having at least 3
digits (base 10).
Please checkout wikipedia for more information.
A circular prime is a prime number with the property that the number
generated at each intermediate step when cyclically permuting its
(base 10) digits will also be prime.
Output
113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933
I use next_prime and is_prime from Math::Prime::Util to produce
primes, and check for primality all (from List::Util) of the cyclic
permutations of succesive prime numbers. I keep a list of seen numbers
to avoid repetitions. This yields a three liner.
perl -MMath::Prime::Util=next_prime,is_prime -MList::Util=all -E '
$n=0; $p=99; while($n<10){ next if $s{$p=next_prime($p)}; @p=split "", $p;
map {$s{$_}=1} @a=(map {join "",@p[$_..@p-1],@p[0..$_-1]} (0..@p-1));
push(@c, $p), ++$n if all {is_prime($_)} @a;} say join " ","Output:", @c'
Results:
Output: 113 197 199 337 1193 3779 11939 19937 193939 199933
The full code is
1 # Perl weekly challenge 167
2 # Task 1: Circular prime
3 #
4 # See https://wlmb.github.io/2022/05/30/PWC167/#task-1-circular-prime
5 use v5.12;
6 use warnings;
7 use Math::Prime::Util qw(next_prime is_prime);
8 use List::Util qw(all);
9 my $count=0;
10 my $candidate=99;
11 my $desired=shift//10; # allow reading number of circular primes from @ARGV
12 say "You are too voracious; next time choose a number <=10" unless $desired <=10;
13 $desired=10 if $desired >10;
14 my %seen;
15 my @found;
16 while($count < $desired){
17 next if $seen{$candidate=next_prime($candidate)};
18 my @digits=split "", $candidate;
19 map {$seen{$_}=1}
20 my @cyclic=map {join "",@digits[$_..@digits-1],@digits[0..$_-1]} (0..@digits-1);
21 push(@found, $candidate), ++$count if all {is_prime($_)} @cyclic;
22 }
23 say join " ","The first $desired decimal circular primes are:", @found;
./ch-1.pl 5
./ch-1.pl 20
Results:
The first 5 decimal circular primes are: 113 197 199 337 1193
You are too optimistic; next time choose a number <=10
The first 10 decimal circular primes are: 113 197 199 337 1193 3779 11939 19937 193939 199933
Task 2: Gamma Function
Submitted by: Mohammad S Anwar
Implement subroutine gamma() using the Lanczos approximation method.
Example
print gamma(3); # 1.99
print gamma(5); # 24
print gamma(7); # 719.99
I could cheat and take advantage of the Gnu Scientific Library implementation of the Lanczos method, as available in the Perl Data Language (PDL) to write a one-liner
Example: Calculate 2!, 4!, and 6!:
perl -MPDL -MPDL::GSLSF::GAMMA -E '
$x=pdl(@ARGV); ($g,$e)=gsl_sf_gamma($x); say "gamma($x)=", $g
' 3 5 7
Results:
gamma([3 5 7])=[2 24 720]
Example: Calculate the gamma functions of a few semiintegers
perl -MPDL -MPDL::GSLSF::GAMMA -E '
$x=pdl(@ARGV); ($g,$e)=gsl_sf_gamma($x); say "gamma($x)=", $g
' .5 1.5 2.5
Results:
gamma([0.5 1.5 2.5])=[ 1.7724539 0.88622693 1.3293404]
The results should have been √π, 0.5×√π, and 1.5×0.5×√π, as verified below,
perl -MPDL -MPDL::GSLSF::GAMMA -E '$sp=sqrt(4*atan2(1,1)); say pdl($sp, .5*$sp, 1.5*.5*$sp);'
Results:
[ 1.7724539 0.88622693 1.3293404]
More challenging though is to actually implement the algorithm. I make a straightforward implementation of the first and fourth equations from the wikipedia article, taking the weights from the page Coefficients for the Lanczos Approximation to the Gamma Function at MROB. I use the Perl Data Language (PDL) to simplify the code avoiding explicit iterations and to allow complex arguments.
1 # Perl weekly challenge 167
2 # Task 2: Gamma function
3 #
4 # See https://wlmb.github.io/2022/05/30/PWC167/#task-2-gamma-function
5 use v5.12;
6 use warnings;
7 use PDL;
8 use PDL::NiceSlice;
9 my $pi=4*atan2(1,1);
10 # coefficients for g=7, n=9, from https://mrob.com/pub/ries/lanczos-gamma.html
11 my $c=pdl[qw( 0.99999999999980993
12 676.5203681218851
13 -1259.1392167224028
14 771.32342877765313
15 -176.61502916214059
16 12.507343278686905
17 -0.13857109526572012
18 9.9843695780195716e-6
19 1.5056327351493116e-7)];
20 my $g=7;
21 die 'Usage: ./ch-2.pl z1 [z2...] to calculate the Gamma function of z1, z2...'
22 unless @ARGV;
23 # Notice that the arguments may be real or complex, of the form 2+3i
24 my $z=pdl join " ", @ARGV;
25 say "gamma$z=", gamma($z);
26
27 sub gamma {
28 my $z = shift;
29 my $gamma=$z->zeroes; # ndarray for results
30 # compute separately for the cases re($z)<.5 and >=.5
31 my ($small_z, $small_gamma)=where($z, $gamma, $z->re <0.5);
32 my ($large_z, $large_gamma)=where($z, $gamma, $z->re >=0.5);
33 gamma_aux(1-$small_z, $small_gamma);
34 $small_gamma.=$pi/(sin($pi*$small_z)*$small_gamma); # reflection formula
35 gamma_aux($large_z, $large_gamma);
36 return $gamma;
37 }
38
39 sub gamma_aux {
40 my ($z, $res)=@_;
41 return if $z->isempty;
42 my $zm1=$z-1;
43 my $den=$zm1->dummy(0)+$c->sequence; # denominator den_{ki}=z_i+c_k
44 $den((0)).=1; # replace den_{0i}=1
45 my $a=($c/$den)->sumover; # sum_k c_k/den_k
46 # compute result and copy into results array
47 $res.=sqrt(2*$pi)*($zm1+$g+1/2)**($zm1+1/2)*exp(-($zm1+$g+1/2))*$a;
48 }
49
Example. Calculate the factorial of a few positive integers
./ch-2.pl 1 2 3 4 5 6 7 8 9 10
The results
gamma[1 2 3 4 5 6 7 8 9]=[1 1 2 6 24 120 720 5040 40320 362880]
are the factorials of 0..9.
Example. Calculate Γ for a few semi-integers
./ch-2.pl -1.5 -.5 .5 1.5
Results:
gamma[-1.5 -0.5 0.5 1.5]=[ 2.3632718 -3.5449077 1.7724539 0.88622693]
We can verify that Γ(0.5)=√π=1.77245385091 and that Γ(z)=(z-1)Γ(z-1), -1.5×2.3632718=-3.5449077, -.5×(-3.5449077)=1.7724539, and .5*1.7724539=0.88622693.
Example. Calculate the factorial of a few complex numbers
./ch-2.pl 1i 1+1i 1-1i 5+3i 5-3i|fmt
Results:
gamma[i 1+i 1-i 5+3i 5-3i]=[-0.15494982830181-0.498015668118356i
0.498015668118356-0.154949828301811i 0.498015668118356+0.154949828301811i
0.0160418827416748-9.43329328975603i 0.0160418827416748+9.43329328975603i]
These results agree with those in the page: Particular values of the gamma function - Wikipedia.