RogerBW's Blog

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.