RabbitFarm
2022-06-05
Circular Primes and Getting Complex
The examples used here are from the weekly challenge problem statement and demonstrate the working solution.
Part 1
Write a script to find out first 10 circular primes having at least 3 digits (base 10).
Solution
use strict;
use warnings;
use boolean;
use Math::Primality qw/is_prime/;
sub is_circular_prime{
my($x, $circular) = @_;
my @digits = split(//, $x);
my @rotations;
for my $i (0 .. @digits - 1){
@digits = (@digits[1 .. @digits - 1], $digits[0]);
my $candidate = join("", @digits) + 0;
push @rotations, $candidate;
return false if !is_prime($candidate);
}
map{$circular->{$_} = -1} @rotations;
return true;
}
sub first_n_circular_primes{
my($n) = @_;
my $i = 100;
my %circular;
my @circular_primes;
{
if(!$circular{$i} && is_circular_prime($i, \%circular)){
push @circular_primes, $i;
}
$i++;
redo if @circular_primes < $n;
}
return @circular_primes;
}
sub first_10_circular_primes{
return first_n_circular_primes(10);
}
MAIN:{
print join(", ", first_10_circular_primes()) . "\n";
}
Sample Run
$ perl perl/ch-1.pl
113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933
Notes
There is a bit of a trick here where we need to disallow repeated use of previous cycles. For example, 199 and 919 and considered to be the same circular prime (we count the first occurrence only) since 919 is a rotation of 199.
I don't ordinarily use a lot of references, especially hash references, in my Perl code
but here it seems appropriate. It makes sense to break the rotating and primality checking
to it's own function but we also need to track all the unique rotations. Wishing to avoid
a global variable, which in this case wouldn't be all that bad anyway, having a single
hash owned by the caller and updated by the primality checking function makes the most
sense to me. The code is arguably cleaner then if we had multiple return values, to
include the rotations. Another option, which would have avoided the use of a reference
and multiple return values would have been to have is_circular_prime
return either
undef
or an array containing the rotations. This would have added a little extra
bookkeeping code to first_n_circular_primes
in order to maintain the master list of all
seen rotations so I considered it, simply as a matter of style, to be just a little less
elegant than the use of the reference.
Part 2
Implement a subroutine gamma() using the Lanczos approximation method.
Solution
use strict;
use warnings;
use POSIX;
use Math::Complex;
use constant EPSILON => 1e-07;
sub lanczos{
my($z) = @_;
my @p = (676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7);
my $y;
$z = new Math::Complex($z);
if(Re($z) < 0.5){
$y = M_PI / (sin(M_PI * $z) * lanczos(1 - $z));
}
else{
$z -= 1;
my $x = 0.99999999999980993;
for my $i (0 .. @p - 1){
$x += ($p[$i] / ($z + $i + 1));
}
my $t = $z + @p - 0.5;
$y = sqrt(2 * M_PI) * $t ** ($z + 0.5) * exp(-1 * $t) * $x;
}
return Re($y) if abs(Im($y)) <= EPSILON;
return $y;
}
sub gamma{
return lanczos(@_);
}
MAIN:{
printf("%.2f\n",gamma(3));
printf("%.2f\n",gamma(5));
printf("%.2f\n",gamma(7));
}
Sample Run
$ perl perl/ch-2.pl
2.00
24.00
720.00
Notes
The code here is based on a Python sample code that accompanies the Wikipedia article and there really isn't much room for additional stylistic flourishes. Well, maybe that for loop could have been a map. For this sort of numeric algorithm there really isn't much variation in what is otherwise a fairly raw computation.
The interesting thing here is that it is by all appearances a faithful representation of
the Lanczos Approximation and yet the answers seem to siffer from a slight floating point
accuracy issue. That is the expected answers vary from what is computed here by a small
decimal part, apparently anyway. Perl is generally quite good at these sorts of things so
getting to the bottom of this may require a bit more investigation! I wonder if it has to
do with how Math::Complex
handles the real part of the number?
References
posted at: 10:46 by: Adam Russell | path: /perl | permanent link to this entry