Merge pull request #6186 from wlmb/challenges

Solve PWC167

3 files changed, 74 insertions, 0 deletions

diff --git a/challenge-167/wlmb/blog.txt b/challenge-167/wlmb/blog.txt
new file mode 100644
index 0000000000..459aeef25b
--- /dev/null
+++ b/challenge-167/wlmb/blog.txt
@@ -0,0 +1 @@
+https://wlmb.github.io/2022/05/30/PWC167/
diff --git a/challenge-167/wlmb/perl/ch-1.pl b/challenge-167/wlmb/perl/ch-1.pl
new file mode 100755
index 0000000000..798af93fee
--- /dev/null
+++ b/challenge-167/wlmb/perl/ch-1.pl
@@ -0,0 +1,24 @@
+#!/usr/bin/env perl
+# Perl weekly challenge 167
+# Task 1: Circular prime
+#
+# See https://wlmb.github.io/2022/05/30/PWC167/#task-1-circular-prime
+use v5.12;
+use warnings;
+use Math::Prime::Util qw(next_prime is_prime);
+use List::Util qw(all);
+my $count=0;
+my $candidate=99;
+my $desired=shift//10; # allow reading number of circular primes from @ARGV
+say "You are too voracious; next time choose a number <=10" unless $desired <=10;
+$desired=10 if $desired >10;
+my %seen;
+my @found;
+while($count < $desired){
+    next if $seen{$candidate=next_prime($candidate)};
+    my @digits=split "", $candidate;
+    map {$seen{$_}=1}
+        my @cyclic=map {join "",@digits[$_..@digits-1],@digits[0..$_-1]} (0..@digits-1);
+    push(@found, $candidate), ++$count if all {is_prime($_)} @cyclic;
+}
+say join " ","The first $desired decimal circular primes are:", @found;
diff --git a/challenge-167/wlmb/perl/ch-2.pl b/challenge-167/wlmb/perl/ch-2.pl
new file mode 100755
index 0000000000..1fb81d137a
--- /dev/null
+++ b/challenge-167/wlmb/perl/ch-2.pl
@@ -0,0 +1,49 @@
+#!/usr/bin/env perl
+# Perl weekly challenge 167
+# Task 2: Gamma function
+#
+# See https://wlmb.github.io/2022/05/30/PWC167/#task-2-gamma-function
+use v5.12;
+use warnings;
+use PDL;
+use PDL::NiceSlice;
+my $pi=4*atan2(1,1);
+# coefficients for g=7, n=9, from https://mrob.com/pub/ries/lanczos-gamma.html
+my $c=pdl[qw(   0.99999999999980993
+              676.5203681218851
+	    -1259.1392167224028
+	       771.32342877765313
+	      -176.61502916214059
+	        12.507343278686905
+	        -0.13857109526572012
+	         9.9843695780195716e-6
+	         1.5056327351493116e-7)];
+my $g=7;
+die 'Usage: ./ch-2.pl z1 [z2...] to calculate the Gamma function of z1, z2...'
+     unless @ARGV;
+# Notice that the arguments may be real or complex, of the form 2+3i
+my $z=pdl join " ", @ARGV;
+say "gamma$z=", gamma($z);
+
+sub gamma {
+    my $z = shift;
+    my $gamma=$z->zeroes; # ndarray for results
+    # compute separately for the cases re($z)<.5 and >=.5
+    my ($small_z, $small_gamma)=where($z, $gamma, $z->re <0.5);
+    my ($large_z, $large_gamma)=where($z, $gamma, $z->re >=0.5);
+    gamma_aux(1-$small_z, $small_gamma);
+    $small_gamma.=$pi/(sin($pi*$small_z)*$small_gamma); # reflection formula
+    gamma_aux($large_z, $large_gamma);
+    return $gamma;
+}
+
+sub gamma_aux {
+    my ($z, $res)=@_;
+    return if $z->isempty;
+    my $zm1=$z-1;
+    my $den=$zm1->dummy(0)+$c->sequence; # denominator den_{ki}=z_i+c_k
+    $den((0)).=1; # replace den_{0i}=1
+    my $a=($c/$den)->sumover; # sum_k  c_k/den_k
+    # compute result and copy into results array
+    $res.=sqrt(2*$pi)*($zm1+$g+1/2)**($zm1+1/2)*exp(-($zm1+$g+1/2))*$a;
+}