# RogerBW's Blog

The Weekly Challenge 158: Prime Cuban Additives 31 March 2022

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

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;
}

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.