Merge pull request #6187 from rjt-pl/master

rjt's week 167 and 138

5 files changed, 140 insertions, 17 deletions

diff --git a/challenge-167/ryan-thompson/README.md b/challenge-167/ryan-thompson/README.md
index dafa7e82f7..5047df481c 100644
--- a/challenge-167/ryan-thompson/README.md
+++ b/challenge-167/ryan-thompson/README.md
@@ -1,28 +1,16 @@
 # Ryan Thompson
 
-## Week 166 Solutions
+## Week 167 Solutions
 
-### Task 1 › Hexadecimal Words
+### Task 1 › Circular Primes
 
  * [Perl](perl/ch-1.pl)
 
-    $ ./ch-1.pl [options]
-
-        --dict=/path/to/dict    Dictionary location
-        --length=8              Target word or phrase length
-        --max-sub=0.2           Ratio of # substitutions / word length
-        --min-length=3          Minimum word length
-        --nopretty              Print hex only, otherwise, pretty print
-
- * [Perl, simplified](perl/ch-1-short.pl)
-
-### Task 2 › K-Directory Diff
+### Task 2 › Lanczos Approximation
 
  * [Perl](perl/ch-2.pl)
 
-    $ ./ch-2.pl dir1 dir2 ...
-
 ## Blogs
 
- * [Hexadecimal Words](https://ry.ca/2022/05/hexadecimal-words/)
- * [K-Directory Diff](https://ry.ca/2022/05/k-directory-diff/)
+ * [Circular Primes](https://ry.ca/2022/05/circular-primes/)
+ * [Lanczos Approximation](https://ry.ca/2022/05/lanczos-approximation/)
diff --git a/challenge-167/ryan-thompson/blog.txt b/challenge-167/ryan-thompson/blog.txt
new file mode 100644
index 0000000000..04da2d0eed
--- /dev/null
+++ b/challenge-167/ryan-thompson/blog.txt
@@ -0,0 +1 @@
+https://ry.ca/2022/05/circular-primes/
diff --git a/challenge-167/ryan-thompson/blog1.txt b/challenge-167/ryan-thompson/blog1.txt
new file mode 100644
index 0000000000..3eff05d9dc
--- /dev/null
+++ b/challenge-167/ryan-thompson/blog1.txt
@@ -0,0 +1 @@
+https://ry.ca/2022/05/lanczos-approximation/
diff --git a/challenge-167/ryan-thompson/perl/ch-1.pl b/challenge-167/ryan-thompson/perl/ch-1.pl
new file mode 100755
index 0000000000..93054d8c75
--- /dev/null
+++ b/challenge-167/ryan-thompson/perl/ch-1.pl
@@ -0,0 +1,49 @@
+#!/usr/bin/env perl
+#
+# ch-1.pl - Circular primes (OEIS A016114)
+#
+# Specifically, the lowest circular prime in each group. So, for example,
+# 113 is a circular prime, as are the rotations 131 and 311, but only 113
+# would be included in the output.
+#
+# 2022 Ryan Thompson <rjt@cpan.org>
+
+use 5.010;
+use warnings;
+use strict;
+no warnings 'uninitialized';
+
+# Just built a sieve a few weeks ago.
+# I have a 3-month Eratosthenes cooldown.
+use Math::Prime::Util qw< primes is_prime >;
+
+my $target = 13; # How many circular primes we want
+my %circ;        # $circ{n} = 1 if n is a circular prime
+
+for (my $digits = 3; ; $digits++) {
+    last if scalar keys %circ >= $target;
+
+    my ($lo, $hi) = (10**($digits-1), 10**$digits-1);
+
+    # All circular primes > 6 digits are repunits
+    if ($digits > 6) {
+        my $n = 1 x $digits;
+        say "R_$digits" if is_prime($n);
+        $circ{$n} = 1;
+        next;
+    }
+
+    my @primes = @{ primes($lo, $hi) };     # All $digits length primes
+    my %prime  = map { $_ => 1 } @primes;
+
+    # Now check every $digits-digit prime for circularity
+N:  for my $n (@primes) {
+        my $rot = $n;
+        for (1..$digits-1) {
+            $rot = chop($rot) . $rot;
+            next N if $circ{$rot} or not $prime{$rot};
+        }
+        say $n;
+        $circ{$n} = 1;
+    }
+}
diff --git a/challenge-167/ryan-thompson/perl/ch-2.pl b/challenge-167/ryan-thompson/perl/ch-2.pl
new file mode 100755
index 0000000000..d0e6558cf9
--- /dev/null
+++ b/challenge-167/ryan-thompson/perl/ch-2.pl
@@ -0,0 +1,84 @@
+#!/usr/bin/env perl
+#
+# ch-2.pl - Gamma function, Lanczos approximation
+#
+# 2022 Ryan Thompson <rjt@cpan.org>
+
+use 5.026;
+use warnings;
+use strict;
+use List::Util qw< sum >;
+use Math::Complex;
+
+use experimental 'signatures';
+no warnings 'uninitialized';
+
+use constant {
+    g   => 7,    # Arbitrary
+   tol  => 1e-7, # Required precision in testing
+};
+
+# These coefficients are only valid for g = 7, @p = 9
+# Credit: https://mrob.com/pub/ries/lanczos-gamma.html
+my @p = (
+     0.99999999999980993227684700473478,
+   676.520368121885098567009190444019,
+ -1259.13921672240287047156078755283,
+   771.3234287776530788486528258894,
+  -176.61502916214059906584551354,
+    12.507343278686904814458936853,
+    -0.13857109526572011689554707,
+     9.984369578019570859563e-6,
+     1.50563273514931155834e-7,
+);
+
+# Lanczos approximation of the gamma function
+# Accepts a real or complex argument
+sub gamma($z) {
+    return pi / (sin($z * pi) * gamma(1 - $z)) if Re($z) < 0.5;
+
+    my $A = $p[0] + sum map { $p[$_] / ($z + $_ - 1) } 1..$#p;
+    my $t = $z + g - 0.5;
+
+    sqrt(2*pi) * $A * $t**($z-0.5) * exp(-$t);
+}
+
+say gamma($ARGV[0]) and exit if defined $ARGV[0];
+
+#
+# Tests
+#
+
+use Test::More;
+
+# Real tests
+my %tests = (
+    # Real tests
+             1 =>  1,
+             2 =>  1,
+             3 =>  2,
+             4 =>  6,
+             5 => 24,
+           0.5 =>       sqrt(pi),
+           1.5 => 0.5 * sqrt(pi),
+          -0.5 =>  -2 * sqrt(pi),
+
+   # Complex tests
+      '1-1i'   => 0.4980156681 + 0.1549498283*i,
+    '0.5+0.5i' => 0.8181639995 - 0.7633138287*i,
+      '5+3i'   => 0.0160418827 - 9.4332932898*i,
+      '5-3i'   => 0.0160418827 + 9.4332932898*i,
+);
+gamma_ok($_, $tests{$_}) for sort keys %tests;
+
+done_testing;
+
+# Test complex or real inputs
+sub gamma_ok($z,$exp) {
+    my $gamma = gamma(cplx($z));
+
+    ok(( abs(Re($gamma) - Re($exp)) <= tol 
+     and abs(Im($gamma) - Im($exp)) <= tol ),
+        "gamma($z) = $exp within " .tol)
+            or diag "  got $gamma";
+}