RogerBW's Blog

The Weekly Challenge 167: Circling the Gamma 02 June 2022

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.

Search
Archive
Tags 1920s 1930s 1940s 1950s 1960s 1970s 1980s 1990s 2000s 2010s 3d printing action aeronautics aikakirja anecdote animation anime army astronomy audio audio tech aviation base commerce battletech beer boardgaming book of the week bookmonth chain of command children chris chronicle church of no redeeming virtues cold war comedy computing contemporary cornish smuggler cosmic encounter coup covid-19 crime cycling dead of winter doctor who documentary drama driving drone ecchi economics espionage essen 2015 essen 2016 essen 2017 essen 2018 essen 2019 essen 2022 existential risk falklands war fandom fanfic fantasy feminism film firefly first world war flash point flight simulation food garmin drive gazebo genesys geocaching geodata gin gkp gurps gurps 101 gus harpoon historical history horror hugo 2014 hugo 2015 hugo 2016 hugo 2017 hugo 2018 hugo 2019 hugo 2020 hugo 2022 hugo-nebula reread in brief avoid instrumented life javascript julian simpson julie enfield kickstarter kotlin learn to play leaving earth linux liquor lovecraftiana lua mecha men with beards museum music mystery naval noir non-fiction one for the brow opera parody paul temple perl perl weekly challenge photography podcast politics postscript powers prediction privacy project woolsack pyracantha python quantum rail raku ranting raspberry pi reading reading boardgames social real life restaurant reviews romance rpg a day rpgs ruby rust science fiction scythe second world war security shipwreck simutrans smartphone south atlantic war squaddies stationery steampunk stuarts suburbia superheroes suspense television the resistance the weekly challenge thirsty meeples thriller tin soldier torg toys trailers travel type 26 type 31 type 45 vietnam war war wargaming weather wives and sweethearts writing about writing x-wing young adult
Special All book reviews, All film reviews
Produced by aikakirja v0.1