I’ve been doing the Weekly
Challenges. The
latest
involved circular primes and the Lanczos approximation. (Note that
this is open until 5 June 2022.)
Task 1: Circular Prime
Write a script to find out first 10 circular primes having at least
3 digits (base 10).
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.
There's a trick here. Yes, 113 is a circular prime - but so are 131
and 311. What's wanted, according to the example, is just the lowest
number in each group.
The Sieve of Eratosthenes relies heavily on having a limit - it will
generate all primes up to that limit with reasonable efficiency, but
it doesn't then have a way of using that list to generate more primes.
So I break this down by digit count: I'll look for all the 3-digit
circular primes, then all the 4-digits, etc., until I have enough.
(Eventually I may develop an extensible version of the sieve code, but
not yet.)
First the cyclic permutation. In all languages, I want to get these
out as numbers, but stringifying the input and manipulating that
seems like the best bet. I realised after I'd written this in all the
languages that it would have been better to double up the string
representation rather than mutating it, then take substrings, thus:
"123" → "123123" → "[123]123", "1[231]23", "12[312]3"
but never mind.
Python:
def cyclicpermute(n):
ss = str(n)
o = []
for p in range(len(ss)):
ss = ss[1:] + ss[0]
o.append(int(ss))
return o
Then the actual search.
def circular(mindigits,ct):
o = []
Initialise the base value with the minimum digits (mindigits = 3 gives
a base of 100).
base = 1
for p in range(mindigits-1):
base = base * 10
while len(o) < ct:
For that number of digits, generate all the primes within it. Put
these in a sorted list, and in a set that I'll use later.
pr = genprimes(base * 10)
prs = set(pr)
Iterate through the list (dropping ones that are lower than the
current base
, which we wastefully calculate multiple times).
for cp in pr:
if cp >= base:
Generate the permutations, and if any of them isn't prime (doesn't
appear in prs
) then drop out and carry on with the next in the list.
v = True
cpp = cyclicpermute(cp)
for cpc in cpp:
if cpc not in prs:
v = False
break
if v:
But if all the permutations are prime: put it on the output list and
exit if we've got enough values. Then mark each of the values in the
permutation as non-prime - so that, having done 113, when we meet 131
we won't also log that one.
o.append(cp)
if len(o) >= ct:
break
for cpc in cpp:
prs.discard(cpc)
If we've done all the n-digit primes without getting enough results,
add a digit and start again.
base = base * 10
return o
Task 2: Gamma Function
Implement subroutine gamma() using the Lanczos approximation method.
Which basically needs complex numbers. I avoid the languages which
don't have easily-accessible complex number libraries (so there's no
JavaScript, Kotlin, Lua or PostScript implementation here). Ruby has
one, but produces very inaccurate results; I'm not sure what's going
on there.
I use the algorithm given (in Python) at
Wikipedia, but
here's the Perl version (using Math::Complex).
Some globals:
my $EPSILON = 1e-7;
my @p = (676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7);
Dropping the imaginary part:
sub drop_imag($z0) {
my $z = $z0;
if (abs(Im($z0)) <= $EPSILON) {
$z = Math::Complex->make(Re($z),0);
}
return $z;
}
The full algorithm using the global coefficients:
sub gamma($z0) {
my $z = Math::Complex->make($z0,0);
my $y;
if (Re($z) < 0.5) {
$y = pi / (sin(pi * $z) * gamma(1-$z));
} else {
$z--;
my $x = 0.99999999999980993;
foreach my $i (0..$#p) {
$x += $p[$i] / ($z + $i + 1);
}
my $t = $z + scalar(@p) - 0.5;
$y = sqrt(2 * pi) * $t ** ($z + 0.5) * exp(-$t) * $x;
}
return drop_imag($y);
}
Full code on
github.
Comments on this post are now closed. If you have particular grounds for adding a late comment, comment on a more recent post quoting the URL of this one.