RogerBW's Blog

I’ve been doing the Weekly Challenges. The latest involved subsets of primes. (Note that this is open until 3 April 2022.)

Task 1: Additive Primes

Write a script to find out all Additive Primes <= 100.

Additive primes are prime numbers for which the sum of their decimal digits are also primes.

I already have code to generate a list of all primes up to a limit, so I'll reuse that. Digit sum is easy enough (repeated modulus, rather than splitting a string which I might do if I were writing only in Perl).

I make use of a relation: the digit sum of a number can be no higher than that number. (I think this is obvious, but I'm not up to producing a formal proof.) So the list of primes that I use to find candidates can be reused as the list of primes that I use to check whether a digit-sum is prime.

In Raku, the new code:

sub digitsum($x0) {
    my $s=0;
    my $x=$x0;
    while ($x > 0) {
        $s += $x % 10;
        $x div= 10;
    }
    return $s;
}

sub additiveprimes($mx) {
    my @o;
    my @p=genprimes($mx);
    my $ps=Set.new(@p);
    for @p -> $q {
        if ($ps{digitsum($q)}:exists) {
            @o.push($q);
        }
    }
    return @o;
}

In some of the languages it was easier to build up the testing-set of primes step by step, during the main loop.

Task 2: First Series Cuban Primes

Write a script to compute first series Cuban Primes <= 1000.

These are the primes that are also a difference between two successive cubes. (E.g. 3³ - 2³ = 27 - 8 = 19.) That can be rearranged into a generating formula: every entry will have the form 3y² + 3y + 1, or rearranged to save an operation, 3 × y × (y+1) + 1. So again I'll use my standard prime generation routine, but this time I just want them in a set, for primality testing of those candidates.

In Python:

def cuban1(mx):
  o=[]

Here's my set of primes.

  ps=set(genprimes(mx))

I don't want to reverse the generation formula to find a stopping point, so in theory I run y all the way up to the limit. Of course the result of the formula is going to be higher than I can use.

  for y in range(1,mx+1):
    q=3*y*(y+1)+1

If the result is higher, we can bail out now. (In fact this will always happen, and the for-loop will never complete.)

    if q > mx:
      break

Failing that, if the number is prime, add it to the output list.

    if q in ps:
      o.append(q)
  return o

Full code on github.